/[escript]/branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
ViewVC logotype

Annotation of /branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3744 - (hide annotations)
Tue Dec 13 06:41:54 2011 UTC (7 years, 10 months ago) by caltinay
File size: 28261 byte(s)
Implemented (Face)Elements->Reduced(Face)Elements interpolation and started
separating DOF and nodes.
Also, implemented operator==() so that a==b for different domain objects a and
b which are in the same state.

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

  ViewVC Help
Powered by ViewVC 1.1.26