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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3754 - (show annotations)
Wed Jan 4 07:07:37 2012 UTC (7 years, 9 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 30637 byte(s)
checkpoint. 3D single PDEs almost working

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 Nodes:
317 case ReducedNodes: //FIXME: reduced
318 case DegreesOfFreedom:
319 case ReducedDegreesOfFreedom: //FIXME: reduced
320 nodesToDOF(target, *const_cast<escript::Data*>(&in));
321 break;
322
323 case Elements:
324 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
325 break;
326
327 case ReducedElements:
328 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
329 break;
330
331 case FaceElements:
332 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
333 break;
334
335 case ReducedFaceElements:
336 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
337 break;
338
339 default:
340 throw RipleyException(msg.str());
341 }
342 break;
343
344 case DegreesOfFreedom:
345 case ReducedDegreesOfFreedom: //FIXME: reduced
346 dofToNodes(target, *const_cast<escript::Data*>(&in));
347 break;
348 default:
349 throw RipleyException(msg.str());
350 }
351 }
352 }
353
354 escript::Data RipleyDomain::getX() const
355 {
356 return escript::continuousFunction(*this).getX();
357 }
358
359 escript::Data RipleyDomain::getNormal() const
360 {
361 return escript::functionOnBoundary(*this).getNormal();
362 }
363
364 escript::Data RipleyDomain::getSize() const
365 {
366 return escript::function(*this).getSize();
367 }
368
369 void RipleyDomain::setToX(escript::Data& arg) const
370 {
371 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
372 *(arg.getFunctionSpace().getDomain()));
373 if (argDomain != *this)
374 throw RipleyException("setToX: Illegal domain of data point locations");
375 if (!arg.isExpanded())
376 throw RipleyException("setToX: Expanded Data object expected");
377
378 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
379 assembleCoordinates(arg);
380 } else {
381 // interpolate the result
382 escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
383 assembleCoordinates(contData);
384 interpolateOnDomain(arg, contData);
385 }
386 }
387
388 bool RipleyDomain::isCellOriented(int fsType) const
389 {
390 switch(fsType) {
391 case Nodes:
392 case DegreesOfFreedom:
393 case ReducedDegreesOfFreedom:
394 return false;
395 case Elements:
396 case FaceElements:
397 case Points:
398 case ReducedElements:
399 case ReducedFaceElements:
400 return true;
401 default:
402 break;
403 }
404 stringstream msg;
405 msg << "isCellOriented(): Illegal function space type " << fsType
406 << " on " << getDescription();
407 throw RipleyException(msg.str());
408 }
409
410 bool RipleyDomain::canTag(int fsType) const
411 {
412 switch(fsType) {
413 case Nodes:
414 case Elements:
415 case ReducedElements:
416 case FaceElements:
417 case ReducedFaceElements:
418 return true;
419 case DegreesOfFreedom:
420 case ReducedDegreesOfFreedom:
421 case Points:
422 return false;
423 default:
424 break;
425 }
426 stringstream msg;
427 msg << "canTag(): Illegal function space type " << fsType << " on "
428 << getDescription();
429 throw RipleyException(msg.str());
430 }
431
432 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
433 {
434 IndexVector* target=NULL;
435 dim_t num=0;
436
437 switch(fsType) {
438 case Nodes:
439 num=getNumNodes();
440 target=&m_nodeTags;
441 break;
442 case Elements:
443 case ReducedElements:
444 num=getNumElements();
445 target=&m_elementTags;
446 break;
447 case FaceElements:
448 case ReducedFaceElements:
449 num=getNumFaceElements();
450 target=&m_faceTags;
451 break;
452 default: {
453 stringstream msg;
454 msg << "setTags(): not implemented for "
455 << functionSpaceTypeAsString(fsType);
456 throw RipleyException(msg.str());
457 }
458 }
459 if (target->size() != num) {
460 target->assign(num, -1);
461 }
462
463 escript::Data& mask=*const_cast<escript::Data*>(&cMask);
464 #pragma omp parallel for
465 for (index_t i=0; i<num; i++) {
466 if (mask.getSampleDataRO(i)[0] > 0) {
467 (*target)[i]=newTag;
468 }
469 }
470 updateTagsInUse(fsType);
471 }
472
473 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
474 {
475 switch(fsType) {
476 case Nodes:
477 if (m_nodeTags.size() > sampleNo)
478 return m_nodeTags[sampleNo];
479 break;
480 case Elements:
481 case ReducedElements:
482 if (m_elementTags.size() > sampleNo)
483 return m_elementTags[sampleNo];
484 break;
485 case FaceElements:
486 case ReducedFaceElements:
487 if (m_faceTags.size() > sampleNo)
488 return m_faceTags[sampleNo];
489 break;
490 default: {
491 stringstream msg;
492 msg << "getTagFromSampleNo(): not implemented for "
493 << functionSpaceTypeAsString(fsType);
494 throw RipleyException(msg.str());
495 }
496 }
497 return -1;
498 }
499
500 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
501 {
502 switch(fsType) {
503 case Nodes:
504 return m_nodeTagsInUse.size();
505 case Elements:
506 case ReducedElements:
507 return m_elementTagsInUse.size();
508 case FaceElements:
509 case ReducedFaceElements:
510 return m_faceTagsInUse.size();
511 default: {
512 stringstream msg;
513 msg << "getNumberOfTagsInUse(): not implemented for "
514 << functionSpaceTypeAsString(fsType);
515 throw RipleyException(msg.str());
516 }
517 }
518 }
519
520 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
521 {
522 switch(fsType) {
523 case Nodes:
524 return &m_nodeTagsInUse[0];
525 case Elements:
526 case ReducedElements:
527 return &m_elementTagsInUse[0];
528 case FaceElements:
529 case ReducedFaceElements:
530 return &m_faceTagsInUse[0];
531 default: {
532 stringstream msg;
533 msg << "borrowListOfTagsInUse(): not implemented for "
534 << functionSpaceTypeAsString(fsType);
535 throw RipleyException(msg.str());
536 }
537 }
538 }
539
540 void RipleyDomain::Print_Mesh_Info(const bool full) const
541 {
542 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
543 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
544 cout << "Number of dimensions: " << m_numDim << endl;
545
546 // write tags
547 if (m_tagMap.size() > 0) {
548 cout << "Tags:" << endl;
549 TagMap::const_iterator it;
550 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
551 cout << " " << setw(5) << it->second << " "
552 << it->first << endl;
553 }
554 }
555 }
556
557 int RipleyDomain::getSystemMatrixTypeId(const int solver,
558 const int preconditioner, const int package, const bool symmetry) const
559 {
560 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
561 package, symmetry, m_mpiInfo);
562 }
563
564 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
565 const int package, const bool symmetry) const
566 {
567 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
568 package, symmetry, m_mpiInfo);
569 }
570
571 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
572 const escript::FunctionSpace& row_functionspace,
573 const int column_blocksize,
574 const escript::FunctionSpace& column_functionspace,
575 const int type) const
576 {
577 bool reduceRowOrder=false;
578 bool reduceColOrder=false;
579 // is the domain right?
580 const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
581 if (row_domain!=*this)
582 throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
583 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
584 if (col_domain!=*this)
585 throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
586 // is the function space type right?
587 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
588 reduceRowOrder=true;
589 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
590 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
591 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
592 reduceColOrder=true;
593 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
594 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
595
596 // generate matrix
597 Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
598 Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
599 row_blocksize, column_blocksize, FALSE);
600 paso::checkPasoError();
601 Paso_SystemMatrixPattern_free(pattern);
602 escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
603 row_functionspace, column_blocksize, column_functionspace));
604 return sma;
605 }
606
607 void RipleyDomain::addPDEToSystem(
608 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
609 const escript::Data& A, const escript::Data& B, const escript::Data& C,
610 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
611 const escript::Data& d, const escript::Data& y,
612 const escript::Data& d_contact, const escript::Data& y_contact,
613 const escript::Data& d_dirac,const escript::Data& y_dirac) const
614 {
615 if (!d_contact.isEmpty() || !y_contact.isEmpty())
616 throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
617
618 paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
619 if (!sma)
620 throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
621
622 if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
623 throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
624
625 //TODO: more input param checks (shape, function space etc)
626
627 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
628
629 if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
630 throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
631
632 const int numEq=S->logical_row_block_size;
633 const int numComp=S->logical_col_block_size;
634
635 if (numEq != numComp)
636 throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
637 //TODO: more system matrix checks
638
639 if (numEq==1)
640 assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
641 else
642 assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
643 }
644
645 void RipleyDomain::setNewX(const escript::Data& arg)
646 {
647 throw RipleyException("setNewX(): Operation not supported");
648 }
649
650 //protected
651 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
652 {
653 const dim_t numComp = in.getDataPointSize();
654 const dim_t dpp = in.getNumDataPointsPerSample();
655 out.requireWrite();
656 #pragma omp parallel for
657 for (index_t i=0; i<in.getNumSamples(); i++) {
658 const double* src = in.getSampleDataRO(i);
659 double* dest = out.getSampleDataRW(i);
660 for (index_t c=0; c<numComp; c++) {
661 double res=0.;
662 for (index_t q=0; q<dpp; q++)
663 res+=src[c+q*numComp];
664 *dest++ = res/dpp;
665 }
666 }
667 }
668
669 //protected
670 void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
671 {
672 const dim_t numComp = in.getDataPointSize();
673 out.requireWrite();
674 #pragma omp parallel for
675 for (index_t i=0; i<in.getNumSamples(); i++) {
676 const double* src = in.getSampleDataRO(i);
677 copy(src, src+numComp, out.getSampleDataRW(i));
678 }
679 }
680
681 //protected
682 void RipleyDomain::updateTagsInUse(int fsType) const
683 {
684 IndexVector* tagsInUse=NULL;
685 const IndexVector* tags=NULL;
686 switch(fsType) {
687 case Nodes:
688 tags=&m_nodeTags;
689 tagsInUse=&m_nodeTagsInUse;
690 break;
691 case Elements:
692 case ReducedElements:
693 tags=&m_elementTags;
694 tagsInUse=&m_elementTagsInUse;
695 break;
696 case FaceElements:
697 case ReducedFaceElements:
698 tags=&m_faceTags;
699 tagsInUse=&m_faceTagsInUse;
700 break;
701 default:
702 return;
703 }
704
705 // gather global unique tag values from tags into tagsInUse
706 tagsInUse->clear();
707 index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
708
709 while (true) {
710 // find smallest value bigger than lastFoundValue
711 minFoundValue = INDEX_T_MAX;
712 #pragma omp parallel private(local_minFoundValue)
713 {
714 local_minFoundValue = minFoundValue;
715 #pragma omp for schedule(static) nowait
716 for (size_t i = 0; i < tags->size(); i++) {
717 const index_t v = (*tags)[i];
718 if ((v > lastFoundValue) && (v < local_minFoundValue))
719 local_minFoundValue = v;
720 }
721 #pragma omp critical
722 {
723 if (local_minFoundValue < minFoundValue)
724 minFoundValue = local_minFoundValue;
725 }
726 }
727 #ifdef ESYS_MPI
728 local_minFoundValue = minFoundValue;
729 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
730 #endif
731
732 // if we found a new value add it to the tagsInUse vector
733 if (minFoundValue < INDEX_T_MAX) {
734 tagsInUse->push_back(minFoundValue);
735 lastFoundValue = minFoundValue;
736 } else
737 break;
738 }
739 }
740
741 //
742 // the methods that follow have to be implemented by the subclasses
743 //
744
745 string RipleyDomain::getDescription() const
746 {
747 throw RipleyException("getDescription() not implemented");
748 }
749
750 void RipleyDomain::write(const string& filename) const
751 {
752 throw RipleyException("write() not implemented");
753 }
754
755 void RipleyDomain::dump(const string& filename) const
756 {
757 throw RipleyException("dump() not implemented");
758 }
759
760 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
761 {
762 throw RipleyException("borrowSampleReferenceIDs() not implemented");
763 }
764
765 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
766 {
767 throw RipleyException("interpolateACross() not implemented");
768 }
769
770 bool RipleyDomain::probeInterpolationACross(int fsType_source,
771 const escript::AbstractDomain&, int fsType_target) const
772 {
773 throw RipleyException("probeInterpolationACross() not implemented");
774 }
775
776 void RipleyDomain::setToNormal(escript::Data& normal) const
777 {
778 throw RipleyException("setToNormal() not implemented");
779 }
780
781 void RipleyDomain::setToSize(escript::Data& size) const
782 {
783 throw RipleyException("setToSize() not implemented");
784 }
785
786 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
787 {
788 throw RipleyException("setToGradient() not implemented");
789 }
790
791 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
792 {
793 throw RipleyException("setToIntegrals() not implemented");
794 }
795
796 bool RipleyDomain::ownSample(int fsType, index_t id) const
797 {
798 throw RipleyException("ownSample() not implemented");
799 }
800
801 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
802 const escript::Data& D, const escript::Data& d,
803 const escript::Data& d_dirac, const bool useHRZ) const
804 {
805 throw RipleyException("addPDEToLumpedSystem() not implemented");
806 }
807
808 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
809 const escript::Data& Y, const escript::Data& y,
810 const escript::Data& y_contact, const escript::Data& y_dirac) const
811 {
812 throw RipleyException("addPDEToRHS() not implemented");
813 }
814
815 void RipleyDomain::addPDEToTransportProblem(
816 escript::AbstractTransportProblem& tp,
817 escript::Data& source, const escript::Data& M,
818 const escript::Data& A, const escript::Data& B, const escript::Data& C,
819 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
820 const escript::Data& d, const escript::Data& y,
821 const escript::Data& d_contact, const escript::Data& y_contact,
822 const escript::Data& d_dirac, const escript::Data& y_dirac) const
823 {
824 throw RipleyException("addPDEToTransportProblem() not implemented");
825 }
826
827 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
828 const int blocksize, const escript::FunctionSpace& functionspace,
829 const int type) const
830 {
831 throw RipleyException("newTransportProblem() not implemented");
832 }
833
834 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
835 bool reducedColOrder) const
836 {
837 throw RipleyException("getPattern() not implemented");
838 }
839
840 dim_t RipleyDomain::getNumDataPointsGlobal() const
841 {
842 throw RipleyException("getNumDataPointsGlobal() not implemented");
843 }
844
845 IndexVector RipleyDomain::getNumNodesPerDim() const
846 {
847 throw RipleyException("getNumNodesPerDim() not implemented");
848 }
849
850 IndexVector RipleyDomain::getNumElementsPerDim() const
851 {
852 throw RipleyException("getNumElementsPerDim() not implemented");
853 }
854
855 IndexVector RipleyDomain::getNumFacesPerBoundary() const
856 {
857 throw RipleyException("getNumFacesPerBoundary() not implemented");
858 }
859
860 IndexVector RipleyDomain::getNodeDistribution() const
861 {
862 throw RipleyException("getNodeDistribution() not implemented");
863 }
864
865 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
866 {
867 throw RipleyException("getFirstCoordAndSpacing() not implemented");
868 }
869
870 dim_t RipleyDomain::getNumFaceElements() const
871 {
872 throw RipleyException("getNumFaceElements() not implemented");
873 }
874
875 dim_t RipleyDomain::getNumElements() const
876 {
877 throw RipleyException("getNumElements() not implemented");
878 }
879
880 dim_t RipleyDomain::getNumNodes() const
881 {
882 throw RipleyException("getNumNodes() not implemented");
883 }
884
885 dim_t RipleyDomain::getNumDOF() const
886 {
887 throw RipleyException("getNumDOF() not implemented");
888 }
889
890 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
891 {
892 throw RipleyException("assembleCoordinates() not implemented");
893 }
894
895 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
896 const escript::Data& A, const escript::Data& B, const escript::Data& C,
897 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
898 const escript::Data& d, const escript::Data& y) const
899 {
900 throw RipleyException("assemblePDESingle() not implemented");
901 }
902
903 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
904 const escript::Data& A, const escript::Data& B, const escript::Data& C,
905 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
906 const escript::Data& d, const escript::Data& y) const
907 {
908 throw RipleyException("assemblePDESystem() not implemented");
909 }
910
911 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
912 {
913 throw RipleyException("interpolateNodesOnElements() not implemented");
914 }
915
916 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
917 {
918 throw RipleyException("interpolateNodesOnFaces() not implemented");
919 }
920
921 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
922 {
923 throw RipleyException("nodesToDOF() not implemented");
924 }
925
926 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
927 {
928 throw RipleyException("dofToNodes() not implemented");
929 }
930
931 } // end of namespace ripley
932

  ViewVC Help
Powered by ViewVC 1.1.26