/[escript]/trunk/ripley/src/RipleyDomain.cpp
ViewVC logotype

Contents of /trunk/ripley/src/RipleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3760 - (show annotations)
Mon Jan 9 05:21:18 2012 UTC (7 years, 8 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 40541 byte(s)
-implemented addPDEToRHS() and setToSize()
-added a few missing calls to requireWrite()
-added assemblePDESystem() to Rectangle but haven't tested yet

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2011 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 #include <ripley/RipleyDomain.h>
15 #include <escript/DataFactory.h>
16 #include <escript/FunctionSpaceFactory.h>
17 #include <pasowrap/SystemMatrixAdapter.h>
18 #include <pasowrap/TransportProblemAdapter.h>
19
20 #include <iomanip>
21
22 using namespace std;
23 using paso::SystemMatrixAdapter;
24 using paso::TransportProblemAdapter;
25
26 namespace ripley {
27
28 escript::Domain_ptr RipleyDomain::loadMesh(const string& filename)
29 {
30 throw RipleyException("loadMesh() not implemented");
31 }
32
33 escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34 {
35 throw RipleyException("readMesh() not implemented");
36 }
37
38 RipleyDomain::RipleyDomain(dim_t dim) :
39 m_numDim(dim),
40 m_status(0)
41 {
42 m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
43 }
44
45 RipleyDomain::~RipleyDomain()
46 {
47 Esys_MPIInfo_free(m_mpiInfo);
48 }
49
50 bool RipleyDomain::operator==(const AbstractDomain& other) const
51 {
52 const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
53 if (o) {
54 return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
55 && m_elementTags==o->m_elementTags
56 && m_faceTags==o->m_faceTags);
57 }
58 return false;
59 }
60
61 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
62 {
63 switch (fsType) {
64 case DegreesOfFreedom:
65 case ReducedDegreesOfFreedom:
66 case Nodes:
67 case ReducedNodes:
68 case Elements:
69 case ReducedElements:
70 case FaceElements:
71 case ReducedFaceElements:
72 case Points:
73 return true;
74 default:
75 break;
76 }
77 return false;
78 }
79
80 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
81 {
82 switch (fsType) {
83 case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84 case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85 case Nodes: return "Ripley_Nodes";
86 case ReducedNodes: return "Ripley_Reduced_Nodes";
87 case Elements: return "Ripley_Elements";
88 case ReducedElements: return "Ripley_Reduced_Elements";
89 case FaceElements: return "Ripley_Face_Elements";
90 case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";
91 case Points: return "Ripley_Points";
92 default:
93 break;
94 }
95 return "Invalid function space type code";
96 }
97
98 pair<int,int> RipleyDomain::getDataShape(int fsType) const
99 {
100 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
101 switch (fsType) {
102 case Nodes:
103 case ReducedNodes: //FIXME: reduced
104 return pair<int,int>(1, getNumNodes());
105 case DegreesOfFreedom:
106 case ReducedDegreesOfFreedom: //FIXME: reduced
107 return pair<int,int>(1, getNumDOF());
108 case Elements:
109 return pair<int,int>(ptsPerSample, getNumElements());
110 case FaceElements:
111 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
112 case ReducedElements:
113 return pair<int,int>(1, getNumElements());
114 case ReducedFaceElements:
115 return pair<int,int>(1, getNumFaceElements());
116 case Points:
117 return pair<int,int>(1, 0); //FIXME: dirac
118 default:
119 break;
120 }
121
122 stringstream msg;
123 msg << "getDataShape(): Unsupported function space type " << fsType
124 << " for " << getDescription();
125 throw RipleyException(msg.str());
126 }
127
128 string RipleyDomain::showTagNames() const
129 {
130 stringstream ret;
131 TagMap::const_iterator it;
132 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
133 if (it!=m_tagMap.begin()) ret << ", ";
134 ret << it->first;
135 }
136 return ret.str();
137 }
138
139 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
140 {
141 /*
142 The idea is to use equivalence classes (i.e. types which can be
143 interpolated back and forth):
144 class 0: DOF <-> Nodes
145 class 1: ReducedDOF <-> ReducedNodes
146 class 2: Points
147 class 3: Elements
148 class 4: ReducedElements
149 class 5: FaceElements
150 class 6: ReducedFaceElements
151
152 There is also a set of lines. Interpolation is possible down a line but not
153 between lines.
154 class 0 and 1 belong to all lines so aren't considered.
155 line 0: class 2
156 line 1: class 3,4
157 line 2: class 5,6
158
159 For classes with multiple members (eg class 1) we have vars to record if
160 there is at least one instance. e.g. hasnodes is true if we have at least
161 one instance of Nodes.
162 */
163 if (fs.empty())
164 return false;
165 vector<bool> hasclass(7, false);
166 vector<int> hasline(3, 0);
167 bool hasnodes=false;
168 bool hasrednodes=false;
169 for (size_t i=0; i<fs.size(); ++i) {
170 switch (fs[i]) {
171 case Nodes: hasnodes=true; // fall through
172 case DegreesOfFreedom:
173 hasclass[0]=true;
174 break;
175 case ReducedNodes: hasrednodes=true; // fall through
176 case ReducedDegreesOfFreedom:
177 hasclass[1]=true;
178 break;
179 case Points:
180 hasline[0]=1;
181 hasclass[2]=true;
182 break;
183 case Elements:
184 hasclass[3]=true;
185 hasline[1]=1;
186 break;
187 case ReducedElements:
188 hasclass[4]=true;
189 hasline[1]=1;
190 break;
191 case FaceElements:
192 hasclass[5]=true;
193 hasline[2]=1;
194 break;
195 case ReducedFaceElements:
196 hasclass[6]=true;
197 hasline[2]=1;
198 break;
199 default:
200 return false;
201 }
202 }
203 int numLines=hasline[0]+hasline[1]+hasline[2];
204
205 // fail if we have more than one leaf group
206 // = there are at least two branches we can't interpolate between
207 if (numLines > 1)
208 return false;
209 else if (numLines==1) {
210 // we have points
211 if (hasline[0]==1)
212 resultcode=Points;
213 else if (hasline[1]==1) {
214 if (hasclass[4])
215 resultcode=ReducedElements;
216 else
217 resultcode=Elements;
218 } else { // hasline[2]==1
219 if (hasclass[6])
220 resultcode=ReducedFaceElements;
221 else
222 resultcode=FaceElements;
223 }
224 } else { // numLines==0
225 if (hasclass[1])
226 // something from class 1
227 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228 else // something from class 0
229 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230 }
231 return true;
232 }
233
234 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
235 int fsType_target) const
236 {
237 if (fsType_target != Nodes &&
238 fsType_target != ReducedNodes &&
239 fsType_target != ReducedDegreesOfFreedom &&
240 fsType_target != DegreesOfFreedom &&
241 fsType_target != Elements &&
242 fsType_target != ReducedElements &&
243 fsType_target != FaceElements &&
244 fsType_target != ReducedFaceElements &&
245 fsType_target != Points) {
246 stringstream msg;
247 msg << "probeInterpolationOnDomain(): Invalid functionspace type "
248 << fsType_target << " for " << getDescription();
249 throw RipleyException(msg.str());
250 }
251
252 switch (fsType_source) {
253 case Nodes:
254 case DegreesOfFreedom:
255 return true;
256 case ReducedNodes:
257 case ReducedDegreesOfFreedom:
258 return (fsType_target != Nodes &&
259 fsType_target != DegreesOfFreedom);
260 case Elements:
261 return (fsType_target==Elements ||
262 fsType_target==ReducedElements);
263 case ReducedElements:
264 return (fsType_target==ReducedElements);
265 case FaceElements:
266 return (fsType_target==FaceElements ||
267 fsType_target==ReducedFaceElements);
268 case ReducedFaceElements:
269 return (fsType_target==ReducedFaceElements);
270 case Points:
271 return (fsType_target==Points);
272
273 default: {
274 stringstream msg;
275 msg << "probeInterpolationOnDomain(): Invalid functionspace type "
276 << fsType_source << " for " << getDescription();
277 throw RipleyException(msg.str());
278 }
279 }
280 }
281
282 void RipleyDomain::interpolateOnDomain(escript::Data& target,
283 const escript::Data& in) const
284 {
285 const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
286 const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
287 if (inDomain != *this)
288 throw RipleyException("Illegal domain of interpolant");
289 if (targetDomain != *this)
290 throw RipleyException("Illegal domain of interpolation target");
291
292 stringstream msg;
293 msg << "interpolateOnDomain() not implemented for function space "
294 << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
295 << " -> "
296 << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
297
298 const int inFS = in.getFunctionSpace().getTypeCode();
299 const int outFS = target.getFunctionSpace().getTypeCode();
300
301 // simplest case: 1:1 copy
302 if (inFS==outFS) {
303 copyData(target, *const_cast<escript::Data*>(&in));
304 // not allowed: reduced->non-reduced
305 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
306 && (outFS==Nodes || outFS==DegreesOfFreedom)) {
307 throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
308 } else if ((inFS==Elements && outFS==ReducedElements)
309 || (inFS==FaceElements && outFS==ReducedFaceElements)) {
310 averageData(target, *const_cast<escript::Data*>(&in));
311 } else {
312 switch (inFS) {
313 case Nodes:
314 case ReducedNodes: //FIXME: reduced
315 switch (outFS) {
316 case DegreesOfFreedom:
317 case ReducedDegreesOfFreedom: //FIXME: reduced
318 if (getMPISize()==1)
319 copyData(target, *const_cast<escript::Data*>(&in));
320 else
321 nodesToDOF(target,*const_cast<escript::Data*>(&in));
322 break;
323
324 case Nodes:
325 case ReducedNodes: //FIXME: reduced
326 copyData(target, *const_cast<escript::Data*>(&in));
327 break;
328
329 case Elements:
330 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
331 break;
332
333 case ReducedElements:
334 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
335 break;
336
337 case FaceElements:
338 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
339 break;
340
341 case ReducedFaceElements:
342 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
343 break;
344
345 default:
346 throw RipleyException(msg.str());
347 }
348 break;
349
350 case DegreesOfFreedom:
351 case ReducedDegreesOfFreedom: //FIXME: reduced
352 switch (outFS) {
353 case Nodes:
354 case ReducedNodes: //FIXME: reduced
355 if (getMPISize()==1)
356 copyData(target, *const_cast<escript::Data*>(&in));
357 else
358 dofToNodes(target, *const_cast<escript::Data*>(&in));
359 break;
360
361 case DegreesOfFreedom:
362 case ReducedDegreesOfFreedom: //FIXME: reduced
363 copyData(target, *const_cast<escript::Data*>(&in));
364 break;
365
366 case Elements:
367 case ReducedElements:
368 if (getMPISize()==1) {
369 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
370 } else {
371 escript::Data contIn=escript::Data(in, continuousFunction(*this));
372 interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
373 }
374 break;
375
376 case FaceElements:
377 case ReducedFaceElements:
378 if (getMPISize()==1) {
379 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
380 } else {
381 escript::Data contIn=escript::Data(in, continuousFunction(*this));
382 interpolateNodesOnElements(target, contIn, outFS==ReducedFaceElements);
383 }
384 break;
385
386 default:
387 throw RipleyException(msg.str());
388 }
389 break;
390 default:
391 throw RipleyException(msg.str());
392 }
393 }
394 }
395
396 escript::Data RipleyDomain::getX() const
397 {
398 return escript::continuousFunction(*this).getX();
399 }
400
401 escript::Data RipleyDomain::getNormal() const
402 {
403 return escript::functionOnBoundary(*this).getNormal();
404 }
405
406 escript::Data RipleyDomain::getSize() const
407 {
408 return escript::function(*this).getSize();
409 }
410
411 void RipleyDomain::setToX(escript::Data& arg) const
412 {
413 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
414 *(arg.getFunctionSpace().getDomain()));
415 if (argDomain != *this)
416 throw RipleyException("setToX: Illegal domain of data point locations");
417 if (!arg.isExpanded())
418 throw RipleyException("setToX: Expanded Data object expected");
419
420 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
421 assembleCoordinates(arg);
422 } else {
423 // interpolate the result
424 escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
425 assembleCoordinates(contData);
426 interpolateOnDomain(arg, contData);
427 }
428 }
429
430 bool RipleyDomain::isCellOriented(int fsType) const
431 {
432 switch(fsType) {
433 case Nodes:
434 case DegreesOfFreedom:
435 case ReducedDegreesOfFreedom:
436 return false;
437 case Elements:
438 case FaceElements:
439 case Points:
440 case ReducedElements:
441 case ReducedFaceElements:
442 return true;
443 default:
444 break;
445 }
446 stringstream msg;
447 msg << "isCellOriented(): Illegal function space type " << fsType
448 << " on " << getDescription();
449 throw RipleyException(msg.str());
450 }
451
452 bool RipleyDomain::canTag(int fsType) const
453 {
454 switch(fsType) {
455 case Nodes:
456 case Elements:
457 case ReducedElements:
458 case FaceElements:
459 case ReducedFaceElements:
460 return true;
461 case DegreesOfFreedom:
462 case ReducedDegreesOfFreedom:
463 case Points:
464 return false;
465 default:
466 break;
467 }
468 stringstream msg;
469 msg << "canTag(): Illegal function space type " << fsType << " on "
470 << getDescription();
471 throw RipleyException(msg.str());
472 }
473
474 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
475 {
476 IndexVector* target=NULL;
477 dim_t num=0;
478
479 switch(fsType) {
480 case Nodes:
481 num=getNumNodes();
482 target=&m_nodeTags;
483 break;
484 case Elements:
485 case ReducedElements:
486 num=getNumElements();
487 target=&m_elementTags;
488 break;
489 case FaceElements:
490 case ReducedFaceElements:
491 num=getNumFaceElements();
492 target=&m_faceTags;
493 break;
494 default: {
495 stringstream msg;
496 msg << "setTags(): not implemented for "
497 << functionSpaceTypeAsString(fsType);
498 throw RipleyException(msg.str());
499 }
500 }
501 if (target->size() != num) {
502 target->assign(num, -1);
503 }
504
505 escript::Data& mask=*const_cast<escript::Data*>(&cMask);
506 #pragma omp parallel for
507 for (index_t i=0; i<num; i++) {
508 if (mask.getSampleDataRO(i)[0] > 0) {
509 (*target)[i]=newTag;
510 }
511 }
512 updateTagsInUse(fsType);
513 }
514
515 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
516 {
517 switch(fsType) {
518 case Nodes:
519 if (m_nodeTags.size() > sampleNo)
520 return m_nodeTags[sampleNo];
521 break;
522 case Elements:
523 case ReducedElements:
524 if (m_elementTags.size() > sampleNo)
525 return m_elementTags[sampleNo];
526 break;
527 case FaceElements:
528 case ReducedFaceElements:
529 if (m_faceTags.size() > sampleNo)
530 return m_faceTags[sampleNo];
531 break;
532 default: {
533 stringstream msg;
534 msg << "getTagFromSampleNo(): not implemented for "
535 << functionSpaceTypeAsString(fsType);
536 throw RipleyException(msg.str());
537 }
538 }
539 return -1;
540 }
541
542 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
543 {
544 switch(fsType) {
545 case Nodes:
546 return m_nodeTagsInUse.size();
547 case Elements:
548 case ReducedElements:
549 return m_elementTagsInUse.size();
550 case FaceElements:
551 case ReducedFaceElements:
552 return m_faceTagsInUse.size();
553 default: {
554 stringstream msg;
555 msg << "getNumberOfTagsInUse(): not implemented for "
556 << functionSpaceTypeAsString(fsType);
557 throw RipleyException(msg.str());
558 }
559 }
560 }
561
562 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
563 {
564 switch(fsType) {
565 case Nodes:
566 return &m_nodeTagsInUse[0];
567 case Elements:
568 case ReducedElements:
569 return &m_elementTagsInUse[0];
570 case FaceElements:
571 case ReducedFaceElements:
572 return &m_faceTagsInUse[0];
573 default: {
574 stringstream msg;
575 msg << "borrowListOfTagsInUse(): not implemented for "
576 << functionSpaceTypeAsString(fsType);
577 throw RipleyException(msg.str());
578 }
579 }
580 }
581
582 void RipleyDomain::Print_Mesh_Info(const bool full) const
583 {
584 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
585 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
586 cout << "Number of dimensions: " << m_numDim << endl;
587
588 // write tags
589 if (m_tagMap.size() > 0) {
590 cout << "Tags:" << endl;
591 TagMap::const_iterator it;
592 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
593 cout << " " << setw(5) << it->second << " "
594 << it->first << endl;
595 }
596 }
597 }
598
599 int RipleyDomain::getSystemMatrixTypeId(const int solver,
600 const int preconditioner, const int package, const bool symmetry) const
601 {
602 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
603 package, symmetry, m_mpiInfo);
604 }
605
606 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
607 const int package, const bool symmetry) const
608 {
609 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
610 package, symmetry, m_mpiInfo);
611 }
612
613 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
614 const escript::FunctionSpace& row_functionspace,
615 const int column_blocksize,
616 const escript::FunctionSpace& column_functionspace,
617 const int type) const
618 {
619 bool reduceRowOrder=false;
620 bool reduceColOrder=false;
621 // is the domain right?
622 const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
623 if (row_domain!=*this)
624 throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
625 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
626 if (col_domain!=*this)
627 throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
628 // is the function space type right?
629 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
630 reduceRowOrder=true;
631 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
632 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
633 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
634 reduceColOrder=true;
635 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
636 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
637
638 // generate matrix
639 Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
640 Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
641 row_blocksize, column_blocksize, FALSE);
642 paso::checkPasoError();
643 Paso_SystemMatrixPattern_free(pattern);
644 escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
645 row_functionspace, column_blocksize, column_functionspace));
646 return sma;
647 }
648
649 void RipleyDomain::addPDEToSystem(
650 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
651 const escript::Data& A, const escript::Data& B, const escript::Data& C,
652 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
653 const escript::Data& d, const escript::Data& y,
654 const escript::Data& d_contact, const escript::Data& y_contact,
655 const escript::Data& d_dirac,const escript::Data& y_dirac) const
656 {
657 if (!d_contact.isEmpty() || !y_contact.isEmpty())
658 throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
659
660 paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
661 if (!sma)
662 throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
663
664 if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
665 throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
666
667 //TODO: more input param checks (shape, function space etc)
668
669 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
670
671 if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
672 throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
673
674 const int numEq=S->logical_row_block_size;
675 const int numComp=S->logical_col_block_size;
676
677 if (numEq != numComp)
678 throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
679 //TODO: more system matrix checks
680
681 if (numEq==1)
682 assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
683 else
684 assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
685 }
686
687 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
688 const escript::Data& Y, const escript::Data& y,
689 const escript::Data& y_contact, const escript::Data& y_dirac) const
690 {
691 if (!y_contact.isEmpty())
692 throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
693
694 if (rhs.isEmpty()) {
695 if (!X.isEmpty() || !Y.isEmpty())
696 throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
697 else
698 return;
699 }
700
701 if (rhs.getDataPointSize() == 1)
702 assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
703 else
704 assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
705 }
706
707 void RipleyDomain::setNewX(const escript::Data& arg)
708 {
709 throw RipleyException("setNewX(): Operation not supported");
710 }
711
712 //protected
713 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
714 {
715 const dim_t numComp = in.getDataPointSize();
716 const dim_t dpp = in.getNumDataPointsPerSample();
717 out.requireWrite();
718 #pragma omp parallel for
719 for (index_t i=0; i<in.getNumSamples(); i++) {
720 const double* src = in.getSampleDataRO(i);
721 double* dest = out.getSampleDataRW(i);
722 for (index_t c=0; c<numComp; c++) {
723 double res=0.;
724 for (index_t q=0; q<dpp; q++)
725 res+=src[c+q*numComp];
726 *dest++ = res/dpp;
727 }
728 }
729 }
730
731 //protected
732 void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
733 {
734 const dim_t numComp = in.getDataPointSize();
735 out.requireWrite();
736 #pragma omp parallel for
737 for (index_t i=0; i<in.getNumSamples(); i++) {
738 const double* src = in.getSampleDataRO(i);
739 copy(src, src+numComp, out.getSampleDataRW(i));
740 }
741 }
742
743 //protected
744 void RipleyDomain::updateTagsInUse(int fsType) const
745 {
746 IndexVector* tagsInUse=NULL;
747 const IndexVector* tags=NULL;
748 switch(fsType) {
749 case Nodes:
750 tags=&m_nodeTags;
751 tagsInUse=&m_nodeTagsInUse;
752 break;
753 case Elements:
754 case ReducedElements:
755 tags=&m_elementTags;
756 tagsInUse=&m_elementTagsInUse;
757 break;
758 case FaceElements:
759 case ReducedFaceElements:
760 tags=&m_faceTags;
761 tagsInUse=&m_faceTagsInUse;
762 break;
763 default:
764 return;
765 }
766
767 // gather global unique tag values from tags into tagsInUse
768 tagsInUse->clear();
769 index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
770
771 while (true) {
772 // find smallest value bigger than lastFoundValue
773 minFoundValue = INDEX_T_MAX;
774 #pragma omp parallel private(local_minFoundValue)
775 {
776 local_minFoundValue = minFoundValue;
777 #pragma omp for schedule(static) nowait
778 for (size_t i = 0; i < tags->size(); i++) {
779 const index_t v = (*tags)[i];
780 if ((v > lastFoundValue) && (v < local_minFoundValue))
781 local_minFoundValue = v;
782 }
783 #pragma omp critical
784 {
785 if (local_minFoundValue < minFoundValue)
786 minFoundValue = local_minFoundValue;
787 }
788 }
789 #ifdef ESYS_MPI
790 local_minFoundValue = minFoundValue;
791 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
792 #endif
793
794 // if we found a new value add it to the tagsInUse vector
795 if (minFoundValue < INDEX_T_MAX) {
796 tagsInUse->push_back(minFoundValue);
797 lastFoundValue = minFoundValue;
798 } else
799 break;
800 }
801 }
802
803 //protected
804 Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
805 const IndexVector& index, const dim_t M, const dim_t N) const
806 {
807 // paso will manage the memory
808 index_t* indexC = MEMALLOC(index.size(), index_t);
809 index_t* ptrC = MEMALLOC(ptr.size(), index_t);
810 copy(index.begin(), index.end(), indexC);
811 copy(ptr.begin(), ptr.end(), ptrC);
812 return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
813 }
814
815 //protected
816 Paso_Pattern* RipleyDomain::createMainPattern() const
817 {
818 IndexVector ptr(1,0);
819 IndexVector index;
820
821 for (index_t i=0; i<getNumDOF(); i++) {
822 // add the DOF itself
823 index.push_back(i);
824 const dim_t num=insertNeighbourNodes(index, i);
825 ptr.push_back(ptr.back()+num+1);
826 }
827
828 return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
829 }
830
831 //protected
832 void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
833 const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
834 {
835 IndexVector ptr(1,0);
836 IndexVector index;
837 for (index_t i=0; i<getNumDOF(); i++) {
838 index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
839 ptr.push_back(ptr.back()+colIndices[i].size());
840 }
841
842 const dim_t M=ptr.size()-1;
843 *colPattern=createPasoPattern(ptr, index, M, N);
844
845 IndexVector rowPtr(1,0);
846 IndexVector rowIndex;
847 for (dim_t id=0; id<N; id++) {
848 dim_t n=0;
849 for (dim_t i=0; i<M; i++) {
850 for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
851 if (index[j]==id) {
852 rowIndex.push_back(i);
853 n++;
854 break;
855 }
856 }
857 }
858 rowPtr.push_back(rowPtr.back()+n);
859 }
860
861 // M and N are now reversed
862 *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
863 }
864
865 //protected
866 void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
867 const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
868 dim_t num_Sol, const vector<double>& array) const
869 {
870 const dim_t numMyCols = mat->pattern->mainPattern->numInput;
871 const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
872 const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
873 const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
874
875 const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
876 const index_t* mainBlock_index = mat->mainBlock->pattern->index;
877 double* mainBlock_val = mat->mainBlock->val;
878 const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
879 const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
880 double* col_coupleBlock_val = mat->col_coupleBlock->val;
881 const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
882 const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
883 double* row_coupleBlock_val = mat->row_coupleBlock->val;
884
885 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
886 // down columns of array
887 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
888 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
889 //printf("row:%d\n", i_row);
890 // only look at the matrix rows stored on this processor
891 if (i_row < numMyRows) {
892 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
893 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
894 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
895 if (i_col < numMyCols) {
896 for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
897 if (mainBlock_index[k] == i_col) {
898 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
899 const dim_t i_Sol=ic+mat->col_block_size*l_col;
900 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
901 const dim_t i_Eq=ir+mat->row_block_size*l_row;
902 mainBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
903 }
904 }
905 break;
906 }
907 }
908 } else {
909 for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
910 //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;
911 if (col_coupleBlock_index[k] == i_col - numMyCols) {
912 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
913 const dim_t i_Sol=ic+mat->col_block_size*l_col;
914 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
915 const dim_t i_Eq=ir+mat->row_block_size*l_row;
916 col_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
917 }
918 }
919 break;
920 }
921 }
922 }
923 }
924 }
925 } else {
926 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
927 // across rows of array
928 for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
929 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
930 if (i_col < numMyCols) {
931 for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
932 k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
933 {
934 //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;
935 if (row_coupleBlock_index[k] == i_col) {
936 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
937 const dim_t i_Sol=ic+mat->col_block_size*l_col;
938 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
939 const dim_t i_Eq=ir+mat->row_block_size*l_row;
940 row_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
941 }
942 }
943 break;
944 }
945 }
946 }
947 }
948 }
949 }
950 }
951 }
952 }
953
954 //
955 // the methods that follow have to be implemented by the subclasses
956 //
957
958 string RipleyDomain::getDescription() const
959 {
960 throw RipleyException("getDescription() not implemented");
961 }
962
963 void RipleyDomain::write(const string& filename) const
964 {
965 throw RipleyException("write() not implemented");
966 }
967
968 void RipleyDomain::dump(const string& filename) const
969 {
970 throw RipleyException("dump() not implemented");
971 }
972
973 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
974 {
975 throw RipleyException("borrowSampleReferenceIDs() not implemented");
976 }
977
978 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
979 {
980 throw RipleyException("interpolateACross() not implemented");
981 }
982
983 bool RipleyDomain::probeInterpolationACross(int fsType_source,
984 const escript::AbstractDomain&, int fsType_target) const
985 {
986 throw RipleyException("probeInterpolationACross() not implemented");
987 }
988
989 void RipleyDomain::setToNormal(escript::Data& normal) const
990 {
991 throw RipleyException("setToNormal() not implemented");
992 }
993
994 void RipleyDomain::setToSize(escript::Data& size) const
995 {
996 throw RipleyException("setToSize() not implemented");
997 }
998
999 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
1000 {
1001 throw RipleyException("setToGradient() not implemented");
1002 }
1003
1004 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
1005 {
1006 throw RipleyException("setToIntegrals() not implemented");
1007 }
1008
1009 bool RipleyDomain::ownSample(int fsType, index_t id) const
1010 {
1011 throw RipleyException("ownSample() not implemented");
1012 }
1013
1014 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1015 const escript::Data& D, const escript::Data& d,
1016 const escript::Data& d_dirac, const bool useHRZ) const
1017 {
1018 throw RipleyException("addPDEToLumpedSystem() not implemented");
1019 }
1020
1021 void RipleyDomain::addPDEToTransportProblem(
1022 escript::AbstractTransportProblem& tp,
1023 escript::Data& source, const escript::Data& M,
1024 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1025 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1026 const escript::Data& d, const escript::Data& y,
1027 const escript::Data& d_contact, const escript::Data& y_contact,
1028 const escript::Data& d_dirac, const escript::Data& y_dirac) const
1029 {
1030 throw RipleyException("addPDEToTransportProblem() not implemented");
1031 }
1032
1033 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
1034 const int blocksize, const escript::FunctionSpace& functionspace,
1035 const int type) const
1036 {
1037 throw RipleyException("newTransportProblem() not implemented");
1038 }
1039
1040 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1041 bool reducedColOrder) const
1042 {
1043 throw RipleyException("getPattern() not implemented");
1044 }
1045
1046 dim_t RipleyDomain::getNumDataPointsGlobal() const
1047 {
1048 throw RipleyException("getNumDataPointsGlobal() not implemented");
1049 }
1050
1051 IndexVector RipleyDomain::getNumNodesPerDim() const
1052 {
1053 throw RipleyException("getNumNodesPerDim() not implemented");
1054 }
1055
1056 IndexVector RipleyDomain::getNumElementsPerDim() const
1057 {
1058 throw RipleyException("getNumElementsPerDim() not implemented");
1059 }
1060
1061 IndexVector RipleyDomain::getNumFacesPerBoundary() const
1062 {
1063 throw RipleyException("getNumFacesPerBoundary() not implemented");
1064 }
1065
1066 IndexVector RipleyDomain::getNodeDistribution() const
1067 {
1068 throw RipleyException("getNodeDistribution() not implemented");
1069 }
1070
1071 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1072 {
1073 throw RipleyException("getFirstCoordAndSpacing() not implemented");
1074 }
1075
1076 dim_t RipleyDomain::getNumFaceElements() const
1077 {
1078 throw RipleyException("getNumFaceElements() not implemented");
1079 }
1080
1081 dim_t RipleyDomain::getNumElements() const
1082 {
1083 throw RipleyException("getNumElements() not implemented");
1084 }
1085
1086 dim_t RipleyDomain::getNumNodes() const
1087 {
1088 throw RipleyException("getNumNodes() not implemented");
1089 }
1090
1091 dim_t RipleyDomain::getNumDOF() const
1092 {
1093 throw RipleyException("getNumDOF() not implemented");
1094 }
1095
1096 dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1097 {
1098 throw RipleyException("insertNeighbourNodes() not implemented");
1099 }
1100
1101 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1102 {
1103 throw RipleyException("assembleCoordinates() not implemented");
1104 }
1105
1106 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1107 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1108 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1109 const escript::Data& d, const escript::Data& y) const
1110 {
1111 throw RipleyException("assemblePDESingle() not implemented");
1112 }
1113
1114 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1115 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1116 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1117 const escript::Data& d, const escript::Data& y) const
1118 {
1119 throw RipleyException("assemblePDESystem() not implemented");
1120 }
1121
1122 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1123 {
1124 throw RipleyException("interpolateNodesOnElements() not implemented");
1125 }
1126
1127 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1128 {
1129 throw RipleyException("interpolateNodesOnFaces() not implemented");
1130 }
1131
1132 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1133 {
1134 throw RipleyException("nodesToDOF() not implemented");
1135 }
1136
1137 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1138 {
1139 throw RipleyException("dofToNodes() not implemented");
1140 }
1141
1142 } // end of namespace ripley
1143

  ViewVC Help
Powered by ViewVC 1.1.26