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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3740 - (hide annotations)
Tue Dec 13 01:37:21 2011 UTC (7 years, 9 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 26502 byte(s)
Merged latest trunk and fixed some smaller issues.

1 caltinay 3670
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 caltinay 3691 #include <escript/FunctionSpaceFactory.h>
17     #include <pasowrap/SystemMatrixAdapter.h>
18     #include <pasowrap/TransportProblemAdapter.h>
19 caltinay 3670
20     #include <iomanip>
21    
22     using namespace std;
23 caltinay 3691 using paso::SystemMatrixAdapter;
24     using paso::TransportProblemAdapter;
25 caltinay 3670
26     namespace ripley {
27    
28 caltinay 3691 escript::Domain_ptr RipleyDomain::loadMesh(const string& filename)
29 caltinay 3670 {
30 caltinay 3691 throw RipleyException("loadMesh() not implemented");
31 caltinay 3670 }
32    
33 caltinay 3691 escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34 caltinay 3670 {
35 caltinay 3691 throw RipleyException("readMesh() not implemented");
36 caltinay 3670 }
37    
38 caltinay 3691 RipleyDomain::RipleyDomain(dim_t dim) :
39     m_numDim(dim),
40     m_status(0)
41 caltinay 3670 {
42 caltinay 3691 m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
43 caltinay 3670 }
44    
45 caltinay 3691 RipleyDomain::~RipleyDomain()
46 caltinay 3670 {
47 caltinay 3691 Esys_MPIInfo_free(m_mpiInfo);
48 caltinay 3670 }
49    
50 caltinay 3691 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
51 caltinay 3670 {
52 caltinay 3691 switch (fsType) {
53 caltinay 3670 case DegreesOfFreedom:
54     case ReducedDegreesOfFreedom:
55     case Nodes:
56     case ReducedNodes:
57     case Elements:
58     case ReducedElements:
59     case FaceElements:
60     case ReducedFaceElements:
61     case Points:
62 caltinay 3691 return true;
63     default:
64 caltinay 3670 break;
65     }
66 caltinay 3691 return false;
67 caltinay 3670 }
68    
69 caltinay 3691 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
70 caltinay 3670 {
71 caltinay 3691 switch (fsType) {
72     case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
73     case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
74     case Nodes: return "Ripley_Nodes";
75     case ReducedNodes: return "Ripley_Reduced_Nodes";
76     case Elements: return "Ripley_Elements";
77     case ReducedElements: return "Ripley_Reduced_Elements";
78     case FaceElements: return "Ripley_Face_Elements";
79     case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";
80     case Points: return "Ripley_Points";
81     default:
82 caltinay 3670 break;
83     }
84 caltinay 3691 return "Invalid function space type code";
85 caltinay 3670 }
86    
87 caltinay 3691 pair<int,int> RipleyDomain::getDataShape(int fsType) const
88 caltinay 3670 {
89 caltinay 3691 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
90     switch (fsType) {
91 caltinay 3670 case Nodes:
92 caltinay 3713 case ReducedNodes:
93 caltinay 3670 case DegreesOfFreedom:
94 caltinay 3722 case ReducedDegreesOfFreedom:
95 caltinay 3691 return pair<int,int>(1, getNumNodes());
96 caltinay 3670 case Elements:
97 caltinay 3691 return pair<int,int>(ptsPerSample, getNumElements());
98 caltinay 3670 case FaceElements:
99 caltinay 3691 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
100 caltinay 3711 case ReducedElements:
101     return pair<int,int>(1, getNumElements());
102     case ReducedFaceElements:
103     return pair<int,int>(1, getNumFaceElements());
104 caltinay 3697 case Points:
105 caltinay 3733 return pair<int,int>(1, 0); //FIXME: dirac
106 caltinay 3670 default:
107     break;
108     }
109    
110 caltinay 3691 stringstream msg;
111 caltinay 3740 msg << "getDataShape(): Unsupported function space type " << fsType
112     << " for " << getDescription();
113 caltinay 3691 throw RipleyException(msg.str());
114 caltinay 3670 }
115    
116 caltinay 3691 string RipleyDomain::showTagNames() const
117 caltinay 3670 {
118 caltinay 3691 stringstream ret;
119     TagMap::const_iterator it;
120     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
121     if (it!=m_tagMap.begin()) ret << ", ";
122     ret << it->first;
123 caltinay 3670 }
124 caltinay 3691 return ret.str();
125 caltinay 3670 }
126    
127     bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
128     {
129 caltinay 3691 /*
130     The idea is to use equivalence classes (i.e. types which can be
131     interpolated back and forth):
132 caltinay 3740 class 0: DOF <-> Nodes
133     class 1: ReducedDOF <-> ReducedNodes
134     class 2: Points
135     class 3: Elements
136     class 4: ReducedElements
137     class 5: FaceElements
138     class 6: ReducedFaceElements
139 caltinay 3670
140 caltinay 3691 There is also a set of lines. Interpolation is possible down a line but not
141     between lines.
142 caltinay 3740 class 0 and 1 belong to all lines so aren't considered.
143     line 0: class 2
144     line 1: class 3,4
145     line 2: class 5,6
146 caltinay 3670
147 caltinay 3740 For classes with multiple members (eg class 1) we have vars to record if
148 caltinay 3670 there is at least one instance. e.g. hasnodes is true if we have at least
149     one instance of Nodes.
150     */
151     if (fs.empty())
152     return false;
153 caltinay 3740 vector<bool> hasclass(7, false);
154     vector<int> hasline(3, 0);
155 caltinay 3670 bool hasnodes=false;
156     bool hasrednodes=false;
157 caltinay 3740 for (size_t i=0; i<fs.size(); ++i) {
158 caltinay 3670 switch (fs[i]) {
159 caltinay 3740 case Nodes: hasnodes=true; // fall through
160 caltinay 3670 case DegreesOfFreedom:
161 caltinay 3740 hasclass[0]=true;
162 caltinay 3670 break;
163 caltinay 3740 case ReducedNodes: hasrednodes=true; // fall through
164 caltinay 3670 case ReducedDegreesOfFreedom:
165 caltinay 3740 hasclass[1]=true;
166 caltinay 3670 break;
167     case Points:
168     hasline[0]=1;
169 caltinay 3740 hasclass[2]=true;
170 caltinay 3670 break;
171     case Elements:
172 caltinay 3740 hasclass[3]=true;
173 caltinay 3670 hasline[1]=1;
174     break;
175     case ReducedElements:
176 caltinay 3740 hasclass[4]=true;
177 caltinay 3670 hasline[1]=1;
178     break;
179     case FaceElements:
180 caltinay 3740 hasclass[5]=true;
181 caltinay 3670 hasline[2]=1;
182     break;
183     case ReducedFaceElements:
184 caltinay 3740 hasclass[6]=true;
185 caltinay 3670 hasline[2]=1;
186     break;
187     default:
188     return false;
189     }
190     }
191 caltinay 3740 int numLines=hasline[0]+hasline[1]+hasline[2];
192 caltinay 3670
193     // fail if we have more than one leaf group
194     // = there are at least two branches we can't interpolate between
195 caltinay 3740 if (numLines > 1)
196 caltinay 3670 return false;
197 caltinay 3740 else if (numLines==1) {
198 caltinay 3670 // we have points
199 caltinay 3691 if (hasline[0]==1)
200 caltinay 3670 resultcode=Points;
201 caltinay 3691 else if (hasline[1]==1) {
202 caltinay 3740 if (hasclass[4])
203 caltinay 3670 resultcode=ReducedElements;
204 caltinay 3691 else
205 caltinay 3670 resultcode=Elements;
206 caltinay 3740 } else { // hasline[2]==1
207     if (hasclass[6])
208 caltinay 3670 resultcode=ReducedFaceElements;
209 caltinay 3691 else
210 caltinay 3670 resultcode=FaceElements;
211 caltinay 3740 }
212     } else { // numLines==0
213     if (hasclass[1])
214     // something from class 1
215 caltinay 3670 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
216 caltinay 3740 else // something from class 0
217 caltinay 3670 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
218     }
219     return true;
220     }
221    
222 caltinay 3732 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
223     int fsType_target) const
224     {
225     if (fsType_target != Nodes &&
226     fsType_target != ReducedNodes &&
227     fsType_target != ReducedDegreesOfFreedom &&
228     fsType_target != DegreesOfFreedom &&
229     fsType_target != Elements &&
230     fsType_target != ReducedElements &&
231     fsType_target != FaceElements &&
232     fsType_target != ReducedFaceElements &&
233     fsType_target != Points) {
234     stringstream msg;
235     msg << "probeInterpolationOnDomain(): Invalid functionspace type "
236 caltinay 3740 << fsType_target << " for " << getDescription();
237 caltinay 3732 throw RipleyException(msg.str());
238     }
239    
240     switch (fsType_source) {
241     case Nodes:
242     case DegreesOfFreedom:
243     return true;
244     case ReducedNodes:
245     case ReducedDegreesOfFreedom:
246     return (fsType_target != Nodes &&
247     fsType_target != DegreesOfFreedom);
248     case Elements:
249     return (fsType_target==Elements ||
250     fsType_target==ReducedElements);
251     case ReducedElements:
252     return (fsType_target==ReducedElements);
253     case FaceElements:
254     return (fsType_target==FaceElements ||
255     fsType_target==ReducedFaceElements);
256     case ReducedFaceElements:
257     return (fsType_target==ReducedFaceElements);
258     case Points:
259     return (fsType_target==Points);
260    
261     default: {
262     stringstream msg;
263     msg << "probeInterpolationOnDomain(): Invalid functionspace type "
264 caltinay 3740 << fsType_source << " for " << getDescription();
265 caltinay 3732 throw RipleyException(msg.str());
266     }
267     }
268     }
269    
270 caltinay 3701 void RipleyDomain::interpolateOnDomain(escript::Data& target,
271     const escript::Data& in) const
272     {
273     const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
274     const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
275     if (inDomain != *this)
276     throw RipleyException("Illegal domain of interpolant");
277     if (targetDomain != *this)
278     throw RipleyException("Illegal domain of interpolation target");
279    
280     stringstream msg;
281     msg << "interpolateOnDomain() not implemented for function space "
282     << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
283     << " -> "
284     << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
285    
286     switch (in.getFunctionSpace().getTypeCode()) {
287     case Nodes:
288 caltinay 3722 case ReducedNodes:
289 caltinay 3701 case DegreesOfFreedom:
290 caltinay 3722 case ReducedDegreesOfFreedom:
291 caltinay 3701 switch (target.getFunctionSpace().getTypeCode()) {
292 caltinay 3722 case Nodes:
293     case ReducedNodes:
294 caltinay 3701 case DegreesOfFreedom:
295 caltinay 3722 case ReducedDegreesOfFreedom:
296     // FIXME
297     copyNodalData(target, *const_cast<escript::Data*>(&in));
298 caltinay 3701 break;
299    
300     case Elements:
301 caltinay 3711 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
302 caltinay 3701 break;
303    
304 caltinay 3711 case ReducedElements:
305     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
306     break;
307    
308 caltinay 3701 case FaceElements:
309 caltinay 3711 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
310 caltinay 3701 break;
311    
312 caltinay 3711 case ReducedFaceElements:
313     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
314     break;
315    
316 caltinay 3701 default:
317     throw RipleyException(msg.str());
318     }
319     break;
320     default:
321     throw RipleyException(msg.str());
322     }
323     }
324    
325 caltinay 3691 escript::Data RipleyDomain::getX() const
326 caltinay 3670 {
327 caltinay 3691 return escript::continuousFunction(*this).getX();
328     }
329    
330     escript::Data RipleyDomain::getNormal() const
331     {
332     return escript::functionOnBoundary(*this).getNormal();
333     }
334    
335     escript::Data RipleyDomain::getSize() const
336     {
337     return escript::function(*this).getSize();
338     }
339    
340     void RipleyDomain::setToX(escript::Data& arg) const
341     {
342     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
343     *(arg.getFunctionSpace().getDomain()));
344     if (argDomain != *this)
345     throw RipleyException("setToX: Illegal domain of data point locations");
346     if (!arg.isExpanded())
347     throw RipleyException("setToX: Expanded Data object expected");
348    
349     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
350     assembleCoordinates(arg);
351     } else {
352     // interpolate the result
353     escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
354     assembleCoordinates(contData);
355     interpolateOnDomain(arg, contData);
356 caltinay 3670 }
357 caltinay 3691 }
358 caltinay 3670
359 caltinay 3691 bool RipleyDomain::isCellOriented(int fsType) const
360     {
361     switch(fsType) {
362 caltinay 3670 case Nodes:
363     case DegreesOfFreedom:
364     case ReducedDegreesOfFreedom:
365 caltinay 3691 return false;
366 caltinay 3670 case Elements:
367 caltinay 3691 case FaceElements:
368     case Points:
369 caltinay 3670 case ReducedElements:
370     case ReducedFaceElements:
371 caltinay 3691 return true;
372 caltinay 3670 default:
373 caltinay 3691 break;
374 caltinay 3670 }
375 caltinay 3691 stringstream msg;
376 caltinay 3740 msg << "isCellOriented(): Illegal function space type " << fsType
377     << " on " << getDescription();
378 caltinay 3691 throw RipleyException(msg.str());
379 caltinay 3670 }
380    
381 caltinay 3722 bool RipleyDomain::canTag(int fsType) const
382     {
383     switch(fsType) {
384     case Nodes:
385     case Elements:
386     case ReducedElements:
387     case FaceElements:
388     case ReducedFaceElements:
389     return true;
390     case DegreesOfFreedom:
391     case ReducedDegreesOfFreedom:
392     case Points:
393     return false;
394     default:
395     break;
396     }
397     stringstream msg;
398 caltinay 3740 msg << "canTag(): Illegal function space type " << fsType << " on "
399     << getDescription();
400 caltinay 3722 throw RipleyException(msg.str());
401     }
402    
403     void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
404     {
405     IndexVector* target=NULL;
406     dim_t num=0;
407    
408     switch(fsType) {
409     case Nodes:
410     num=getNumNodes();
411     target=&m_nodeTags;
412     break;
413     case Elements:
414     case ReducedElements:
415     num=getNumElements();
416     target=&m_elementTags;
417     break;
418     case FaceElements:
419     case ReducedFaceElements:
420     num=getNumFaceElements();
421     target=&m_faceTags;
422     break;
423     default: {
424     stringstream msg;
425 caltinay 3740 msg << "setTags(): not implemented for "
426     << functionSpaceTypeAsString(fsType);
427 caltinay 3722 throw RipleyException(msg.str());
428     }
429     }
430     if (target->size() != num) {
431     target->assign(num, -1);
432     }
433    
434     escript::Data& mask=*const_cast<escript::Data*>(&cMask);
435     #pragma omp parallel for
436     for (index_t i=0; i<num; i++) {
437     if (mask.getSampleDataRO(i)[0] > 0) {
438     (*target)[i]=newTag;
439     }
440     }
441     updateTagsInUse(fsType);
442     }
443    
444     int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
445     {
446     switch(fsType) {
447     case Nodes:
448     if (m_nodeTags.size() > sampleNo)
449     return m_nodeTags[sampleNo];
450     break;
451     case Elements:
452     case ReducedElements:
453     if (m_elementTags.size() > sampleNo)
454     return m_elementTags[sampleNo];
455     break;
456     case FaceElements:
457     case ReducedFaceElements:
458     if (m_faceTags.size() > sampleNo)
459     return m_faceTags[sampleNo];
460     break;
461     default: {
462     stringstream msg;
463 caltinay 3740 msg << "getTagFromSampleNo(): not implemented for "
464     << functionSpaceTypeAsString(fsType);
465 caltinay 3722 throw RipleyException(msg.str());
466     }
467     }
468     return -1;
469     }
470    
471     int RipleyDomain::getNumberOfTagsInUse(int fsType) const
472     {
473     switch(fsType) {
474     case Nodes:
475     return m_nodeTagsInUse.size();
476     case Elements:
477     case ReducedElements:
478     return m_elementTagsInUse.size();
479     case FaceElements:
480     case ReducedFaceElements:
481     return m_faceTagsInUse.size();
482     default: {
483     stringstream msg;
484 caltinay 3740 msg << "getNumberOfTagsInUse(): not implemented for "
485     << functionSpaceTypeAsString(fsType);
486 caltinay 3722 throw RipleyException(msg.str());
487     }
488     }
489     }
490    
491     const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
492     {
493     switch(fsType) {
494     case Nodes:
495     return &m_nodeTagsInUse[0];
496     case Elements:
497     case ReducedElements:
498     return &m_elementTagsInUse[0];
499     case FaceElements:
500     case ReducedFaceElements:
501     return &m_faceTagsInUse[0];
502     default: {
503     stringstream msg;
504 caltinay 3740 msg << "borrowListOfTagsInUse(): not implemented for "
505     << functionSpaceTypeAsString(fsType);
506 caltinay 3722 throw RipleyException(msg.str());
507     }
508     }
509     }
510    
511 caltinay 3691 void RipleyDomain::Print_Mesh_Info(const bool full) const
512 caltinay 3670 {
513 caltinay 3691 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
514     << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
515     cout << "Number of dimensions: " << m_numDim << endl;
516 caltinay 3670
517 caltinay 3691 // write tags
518     if (m_tagMap.size() > 0) {
519     cout << "Tags:" << endl;
520     TagMap::const_iterator it;
521     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
522     cout << " " << setw(5) << it->second << " "
523     << it->first << endl;
524     }
525 caltinay 3670 }
526     }
527    
528     int RipleyDomain::getSystemMatrixTypeId(const int solver,
529     const int preconditioner, const int package, const bool symmetry) const
530     {
531 caltinay 3691 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
532     package, symmetry, m_mpiInfo);
533 caltinay 3670 }
534    
535     int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
536     const int package, const bool symmetry) const
537     {
538 caltinay 3691 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
539     package, symmetry, m_mpiInfo);
540 caltinay 3670 }
541    
542 caltinay 3691 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
543     const escript::FunctionSpace& row_functionspace,
544     const int column_blocksize,
545     const escript::FunctionSpace& column_functionspace,
546     const int type) const
547 caltinay 3670 {
548 caltinay 3691 bool reduceRowOrder=false;
549     bool reduceColOrder=false;
550     // is the domain right?
551     const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
552     if (row_domain!=*this)
553 caltinay 3740 throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
554 caltinay 3691 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
555     if (col_domain!=*this)
556 caltinay 3740 throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
557 caltinay 3691 // is the function space type right?
558     if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
559     reduceRowOrder=true;
560     else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
561 caltinay 3740 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
562 caltinay 3691 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
563     reduceColOrder=true;
564     else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
565 caltinay 3740 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
566 caltinay 3691
567     // generate matrix
568     Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
569     Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
570     row_blocksize, column_blocksize, FALSE);
571     paso::checkPasoError();
572     Paso_SystemMatrixPattern_free(pattern);
573     escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
574     row_functionspace, column_blocksize, column_functionspace));
575     return sma;
576 caltinay 3670 }
577    
578 caltinay 3697 void RipleyDomain::setNewX(const escript::Data& arg)
579     {
580     throw RipleyException("setNewX(): Operation not supported");
581     }
582    
583 caltinay 3722 //protected
584     void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const
585     {
586     const dim_t numComp = in.getDataPointSize();
587 caltinay 3740 out.requireWrite();
588 caltinay 3722 #pragma omp parallel for
589     for (index_t i=0; i<in.getNumSamples(); i++) {
590     const double* src = in.getSampleDataRO(i);
591     copy(src, src+numComp, out.getSampleDataRW(i));
592     }
593     }
594    
595     //protected
596     void RipleyDomain::updateTagsInUse(int fsType) const
597     {
598     IndexVector* tagsInUse=NULL;
599     const IndexVector* tags=NULL;
600     switch(fsType) {
601     case Nodes:
602     tags=&m_nodeTags;
603     tagsInUse=&m_nodeTagsInUse;
604     break;
605     case Elements:
606     case ReducedElements:
607     tags=&m_elementTags;
608     tagsInUse=&m_elementTagsInUse;
609     break;
610     case FaceElements:
611     case ReducedFaceElements:
612     tags=&m_faceTags;
613     tagsInUse=&m_faceTagsInUse;
614     break;
615     default:
616     return;
617     }
618    
619     // gather global unique tag values from tags into tagsInUse
620     tagsInUse->clear();
621     index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
622    
623     while (true) {
624     // find smallest value bigger than lastFoundValue
625     minFoundValue = INDEX_T_MAX;
626     #pragma omp parallel private(local_minFoundValue)
627     {
628     local_minFoundValue = minFoundValue;
629     #pragma omp for schedule(static) nowait
630     for (size_t i = 0; i < tags->size(); i++) {
631     const index_t v = (*tags)[i];
632     if ((v > lastFoundValue) && (v < local_minFoundValue))
633     local_minFoundValue = v;
634     }
635     #pragma omp critical
636     {
637     if (local_minFoundValue < minFoundValue)
638     minFoundValue = local_minFoundValue;
639     }
640     }
641     #ifdef ESYS_MPI
642     local_minFoundValue = minFoundValue;
643     MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
644     #endif
645 caltinay 3740
646 caltinay 3722 // if we found a new value add it to the tagsInUse vector
647     if (minFoundValue < INDEX_T_MAX) {
648     tagsInUse->push_back(minFoundValue);
649     lastFoundValue = minFoundValue;
650     } else
651     break;
652     }
653     }
654    
655 caltinay 3691 //
656     // the methods that follow have to be implemented by the subclasses
657     //
658    
659     string RipleyDomain::getDescription() const
660 caltinay 3670 {
661 caltinay 3691 throw RipleyException("getDescription() not implemented");
662 caltinay 3670 }
663    
664 caltinay 3691 bool RipleyDomain::operator==(const AbstractDomain& other) const
665 caltinay 3670 {
666 caltinay 3691 throw RipleyException("operator==() not implemented");
667 caltinay 3670 }
668    
669 caltinay 3691 void RipleyDomain::write(const string& filename) const
670 caltinay 3670 {
671 caltinay 3691 throw RipleyException("write() not implemented");
672 caltinay 3670 }
673    
674 caltinay 3691 void RipleyDomain::dump(const string& filename) const
675 caltinay 3670 {
676 caltinay 3691 throw RipleyException("dump() not implemented");
677 caltinay 3670 }
678    
679 caltinay 3691 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
680     {
681     throw RipleyException("borrowSampleReferenceIDs() not implemented");
682     }
683 caltinay 3670
684 caltinay 3691 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
685 caltinay 3670 {
686 caltinay 3691 throw RipleyException("interpolateACross() not implemented");
687 caltinay 3670 }
688    
689 caltinay 3691 bool RipleyDomain::probeInterpolationACross(int fsType_source,
690     const escript::AbstractDomain&, int fsType_target) const
691 caltinay 3670 {
692 caltinay 3691 throw RipleyException("probeInterpolationACross() not implemented");
693 caltinay 3670 }
694    
695 caltinay 3691 void RipleyDomain::setToNormal(escript::Data& normal) const
696 caltinay 3670 {
697 caltinay 3691 throw RipleyException("setToNormal() not implemented");
698 caltinay 3670 }
699    
700 caltinay 3691 void RipleyDomain::setToSize(escript::Data& size) const
701 caltinay 3670 {
702 caltinay 3691 throw RipleyException("setToSize() not implemented");
703 caltinay 3670 }
704    
705 caltinay 3691 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
706 caltinay 3670 {
707 caltinay 3691 throw RipleyException("setToGradient() not implemented");
708 caltinay 3670 }
709    
710 caltinay 3691 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
711 caltinay 3670 {
712 caltinay 3691 throw RipleyException("setToIntegrals() not implemented");
713 caltinay 3670 }
714    
715 caltinay 3691 bool RipleyDomain::ownSample(int fsType, index_t id) const
716 caltinay 3670 {
717 caltinay 3691 throw RipleyException("ownSample() not implemented");
718 caltinay 3670 }
719    
720 caltinay 3691 void RipleyDomain::addPDEToSystem(
721     escript::AbstractSystemMatrix& mat, escript::Data& rhs,
722     const escript::Data& A, const escript::Data& B, const escript::Data& C,
723     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
724     const escript::Data& d, const escript::Data& y,
725     const escript::Data& d_contact, const escript::Data& y_contact,
726     const escript::Data& d_dirac,const escript::Data& y_dirac) const
727     {
728     throw RipleyException("addPDEToSystem() not implemented");
729     }
730 caltinay 3670
731 caltinay 3691 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
732     const escript::Data& D, const escript::Data& d,
733     const escript::Data& d_dirac, const bool useHRZ) const
734     {
735     throw RipleyException("addPDEToLumpedSystem() not implemented");
736 caltinay 3670 }
737    
738 caltinay 3691 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
739     const escript::Data& Y, const escript::Data& y,
740     const escript::Data& y_contact, const escript::Data& y_dirac) const
741 caltinay 3670 {
742 caltinay 3691 throw RipleyException("addPDEToRHS() not implemented");
743 caltinay 3670 }
744    
745 caltinay 3691 void RipleyDomain::addPDEToTransportProblem(
746     escript::AbstractTransportProblem& tp,
747     escript::Data& source, const escript::Data& M,
748     const escript::Data& A, const escript::Data& B, const escript::Data& C,
749     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
750     const escript::Data& d, const escript::Data& y,
751     const escript::Data& d_contact, const escript::Data& y_contact,
752     const escript::Data& d_dirac, const escript::Data& y_dirac) const
753 caltinay 3670 {
754 caltinay 3691 throw RipleyException("addPDEToTransportProblem() not implemented");
755 caltinay 3670 }
756    
757 caltinay 3691 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
758     const int blocksize, const escript::FunctionSpace& functionspace,
759     const int type) const
760 caltinay 3670 {
761 caltinay 3691 throw RipleyException("newTransportProblem() not implemented");
762 caltinay 3670 }
763    
764 caltinay 3691 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
765     bool reducedColOrder) const
766 caltinay 3670 {
767 caltinay 3691 throw RipleyException("getPattern() not implemented");
768 caltinay 3670 }
769    
770 caltinay 3691 dim_t RipleyDomain::getNumDataPointsGlobal() const
771 caltinay 3670 {
772 caltinay 3691 throw RipleyException("getNumDataPointsGlobal() not implemented");
773 caltinay 3670 }
774    
775 caltinay 3697 IndexVector RipleyDomain::getNumNodesPerDim() const
776     {
777     throw RipleyException("getNumNodesPerDim() not implemented");
778     }
779    
780     IndexVector RipleyDomain::getNumElementsPerDim() const
781     {
782     throw RipleyException("getNumElementsPerDim() not implemented");
783     }
784    
785     IndexVector RipleyDomain::getNumFacesPerBoundary() const
786     {
787     throw RipleyException("getNumFacesPerBoundary() not implemented");
788     }
789    
790     IndexVector RipleyDomain::getNodeDistribution() const
791     {
792     throw RipleyException("getNodeDistribution() not implemented");
793     }
794    
795     pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
796     {
797     throw RipleyException("getFirstCoordAndSpacing() not implemented");
798     }
799    
800 caltinay 3691 dim_t RipleyDomain::getNumFaceElements() const
801 caltinay 3670 {
802 caltinay 3691 throw RipleyException("getNumFaceElements() not implemented");
803 caltinay 3670 }
804    
805 caltinay 3691 dim_t RipleyDomain::getNumElements() const
806 caltinay 3670 {
807 caltinay 3691 throw RipleyException("getNumElements() not implemented");
808 caltinay 3670 }
809    
810 caltinay 3691 dim_t RipleyDomain::getNumNodes() const
811 caltinay 3670 {
812 caltinay 3691 throw RipleyException("getNumNodes() not implemented");
813 caltinay 3670 }
814    
815 caltinay 3691 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
816 caltinay 3670 {
817 caltinay 3691 throw RipleyException("assembleCoordinates() not implemented");
818 caltinay 3670 }
819    
820 caltinay 3711 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
821 caltinay 3701 {
822     throw RipleyException("interpolateNodesOnElements() not implemented");
823     }
824 caltinay 3691
825 caltinay 3711 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
826 caltinay 3701 {
827     throw RipleyException("interpolateNodesOnFaces() not implemented");
828     }
829    
830    
831 caltinay 3670 } // end of namespace ripley
832    

  ViewVC Help
Powered by ViewVC 1.1.26