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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3756 - (hide annotations)
Fri Jan 6 02:35:19 2012 UTC (7 years, 8 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 37674 byte(s)
Fixed interpolation from DOF to nodes and moved common code to the base class.

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 3750 case DegreesOfFreedom:
65     case ReducedDegreesOfFreedom:
66 caltinay 3670 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 caltinay 3750 case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84     case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85 caltinay 3691 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 3750 case DegreesOfFreedom:
106     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 3750 class 0: DOF <-> Nodes
145     class 1: ReducedDOF <-> ReducedNodes
146 caltinay 3740 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 3750
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 caltinay 3670 */
163     if (fs.empty())
164     return false;
165 caltinay 3740 vector<bool> hasclass(7, false);
166     vector<int> hasline(3, 0);
167 caltinay 3750 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 3750 case Nodes: hasnodes=true; // fall through
172     case DegreesOfFreedom:
173 caltinay 3740 hasclass[0]=true;
174 caltinay 3670 break;
175 caltinay 3750 case ReducedNodes: hasrednodes=true; // fall through
176     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 3750 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 3750 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 3750 // 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 caltinay 3750 // something from class 1
227     resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228     else // something from class 0
229     resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230 caltinay 3670 }
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 caltinay 3750 fsType_target != ReducedDegreesOfFreedom &&
240     fsType_target != DegreesOfFreedom &&
241 caltinay 3732 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 caltinay 3750 case DegreesOfFreedom:
255 caltinay 3732 return true;
256     case ReducedNodes:
257 caltinay 3750 case ReducedDegreesOfFreedom:
258     return (fsType_target != Nodes &&
259     fsType_target != DegreesOfFreedom);
260 caltinay 3732 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 caltinay 3750 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
306     && (outFS==Nodes || outFS==DegreesOfFreedom)) {
307 caltinay 3744 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 caltinay 3750 case ReducedNodes: //FIXME: reduced
315 caltinay 3744 switch (outFS) {
316 caltinay 3750 case DegreesOfFreedom:
317     case ReducedDegreesOfFreedom: //FIXME: reduced
318 caltinay 3755 if (getMPISize()==1)
319     copyData(target, *const_cast<escript::Data*>(&in));
320     else
321     nodesToDOF(target,*const_cast<escript::Data*>(&in));
322 caltinay 3744 break;
323 caltinay 3701
324 caltinay 3755 case Nodes:
325     case ReducedNodes: //FIXME: reduced
326     copyData(target, *const_cast<escript::Data*>(&in));
327     break;
328    
329 caltinay 3744 case Elements:
330     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
331     break;
332 caltinay 3711
333 caltinay 3744 case ReducedElements:
334     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
335     break;
336 caltinay 3701
337 caltinay 3744 case FaceElements:
338     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
339     break;
340 caltinay 3711
341 caltinay 3744 case ReducedFaceElements:
342     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
343     break;
344    
345     default:
346     throw RipleyException(msg.str());
347     }
348     break;
349 caltinay 3750
350     case DegreesOfFreedom:
351     case ReducedDegreesOfFreedom: //FIXME: reduced
352 caltinay 3755 switch (outFS) {
353     case Nodes:
354     case ReducedNodes: //FIXME: reduced
355     if (getMPISize()==1)
356     copyData(target, *const_cast<escript::Data*>(&in));
357     else
358     dofToNodes(target, *const_cast<escript::Data*>(&in));
359     break;
360    
361     case DegreesOfFreedom:
362     case ReducedDegreesOfFreedom: //FIXME: reduced
363     copyData(target, *const_cast<escript::Data*>(&in));
364     break;
365    
366     case Elements:
367     case ReducedElements:
368     if (getMPISize()==1) {
369     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
370     } else {
371     escript::Data contIn=escript::Data(in, continuousFunction(*this));
372     interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
373     }
374     break;
375    
376     case FaceElements:
377     case ReducedFaceElements:
378     if (getMPISize()==1) {
379     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
380     } else {
381     escript::Data contIn=escript::Data(in, continuousFunction(*this));
382     interpolateNodesOnElements(target, contIn, outFS==ReducedFaceElements);
383     }
384     break;
385    
386     default:
387     throw RipleyException(msg.str());
388     }
389 caltinay 3754 break;
390 caltinay 3744 default:
391     throw RipleyException(msg.str());
392     }
393 caltinay 3701 }
394     }
395    
396 caltinay 3691 escript::Data RipleyDomain::getX() const
397 caltinay 3670 {
398 caltinay 3691 return escript::continuousFunction(*this).getX();
399     }
400    
401     escript::Data RipleyDomain::getNormal() const
402     {
403     return escript::functionOnBoundary(*this).getNormal();
404     }
405    
406     escript::Data RipleyDomain::getSize() const
407     {
408     return escript::function(*this).getSize();
409     }
410    
411     void RipleyDomain::setToX(escript::Data& arg) const
412     {
413     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
414     *(arg.getFunctionSpace().getDomain()));
415     if (argDomain != *this)
416     throw RipleyException("setToX: Illegal domain of data point locations");
417     if (!arg.isExpanded())
418     throw RipleyException("setToX: Expanded Data object expected");
419    
420     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
421     assembleCoordinates(arg);
422     } else {
423     // interpolate the result
424     escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
425     assembleCoordinates(contData);
426     interpolateOnDomain(arg, contData);
427 caltinay 3670 }
428 caltinay 3691 }
429 caltinay 3670
430 caltinay 3691 bool RipleyDomain::isCellOriented(int fsType) const
431     {
432     switch(fsType) {
433 caltinay 3670 case Nodes:
434 caltinay 3750 case DegreesOfFreedom:
435     case ReducedDegreesOfFreedom:
436 caltinay 3691 return false;
437 caltinay 3670 case Elements:
438 caltinay 3691 case FaceElements:
439     case Points:
440 caltinay 3670 case ReducedElements:
441     case ReducedFaceElements:
442 caltinay 3691 return true;
443 caltinay 3670 default:
444 caltinay 3691 break;
445 caltinay 3670 }
446 caltinay 3691 stringstream msg;
447 caltinay 3740 msg << "isCellOriented(): Illegal function space type " << fsType
448     << " on " << getDescription();
449 caltinay 3691 throw RipleyException(msg.str());
450 caltinay 3670 }
451    
452 caltinay 3722 bool RipleyDomain::canTag(int fsType) const
453     {
454     switch(fsType) {
455     case Nodes:
456     case Elements:
457     case ReducedElements:
458     case FaceElements:
459     case ReducedFaceElements:
460     return true;
461 caltinay 3750 case DegreesOfFreedom:
462     case ReducedDegreesOfFreedom:
463 caltinay 3722 case Points:
464     return false;
465     default:
466     break;
467     }
468     stringstream msg;
469 caltinay 3740 msg << "canTag(): Illegal function space type " << fsType << " on "
470     << getDescription();
471 caltinay 3722 throw RipleyException(msg.str());
472     }
473    
474     void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
475     {
476     IndexVector* target=NULL;
477     dim_t num=0;
478    
479     switch(fsType) {
480     case Nodes:
481     num=getNumNodes();
482     target=&m_nodeTags;
483     break;
484     case Elements:
485     case ReducedElements:
486     num=getNumElements();
487     target=&m_elementTags;
488     break;
489     case FaceElements:
490     case ReducedFaceElements:
491     num=getNumFaceElements();
492     target=&m_faceTags;
493     break;
494     default: {
495     stringstream msg;
496 caltinay 3740 msg << "setTags(): not implemented for "
497     << functionSpaceTypeAsString(fsType);
498 caltinay 3722 throw RipleyException(msg.str());
499     }
500     }
501     if (target->size() != num) {
502     target->assign(num, -1);
503     }
504    
505     escript::Data& mask=*const_cast<escript::Data*>(&cMask);
506     #pragma omp parallel for
507     for (index_t i=0; i<num; i++) {
508     if (mask.getSampleDataRO(i)[0] > 0) {
509     (*target)[i]=newTag;
510     }
511     }
512     updateTagsInUse(fsType);
513     }
514    
515     int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
516     {
517     switch(fsType) {
518     case Nodes:
519     if (m_nodeTags.size() > sampleNo)
520     return m_nodeTags[sampleNo];
521     break;
522     case Elements:
523     case ReducedElements:
524     if (m_elementTags.size() > sampleNo)
525     return m_elementTags[sampleNo];
526     break;
527     case FaceElements:
528     case ReducedFaceElements:
529     if (m_faceTags.size() > sampleNo)
530     return m_faceTags[sampleNo];
531     break;
532     default: {
533     stringstream msg;
534 caltinay 3740 msg << "getTagFromSampleNo(): not implemented for "
535     << functionSpaceTypeAsString(fsType);
536 caltinay 3722 throw RipleyException(msg.str());
537     }
538     }
539     return -1;
540     }
541    
542     int RipleyDomain::getNumberOfTagsInUse(int fsType) const
543     {
544     switch(fsType) {
545     case Nodes:
546     return m_nodeTagsInUse.size();
547     case Elements:
548     case ReducedElements:
549     return m_elementTagsInUse.size();
550     case FaceElements:
551     case ReducedFaceElements:
552     return m_faceTagsInUse.size();
553     default: {
554     stringstream msg;
555 caltinay 3740 msg << "getNumberOfTagsInUse(): not implemented for "
556     << functionSpaceTypeAsString(fsType);
557 caltinay 3722 throw RipleyException(msg.str());
558     }
559     }
560     }
561    
562     const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
563     {
564     switch(fsType) {
565     case Nodes:
566     return &m_nodeTagsInUse[0];
567     case Elements:
568     case ReducedElements:
569     return &m_elementTagsInUse[0];
570     case FaceElements:
571     case ReducedFaceElements:
572     return &m_faceTagsInUse[0];
573     default: {
574     stringstream msg;
575 caltinay 3740 msg << "borrowListOfTagsInUse(): not implemented for "
576     << functionSpaceTypeAsString(fsType);
577 caltinay 3722 throw RipleyException(msg.str());
578     }
579     }
580     }
581    
582 caltinay 3691 void RipleyDomain::Print_Mesh_Info(const bool full) const
583 caltinay 3670 {
584 caltinay 3691 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
585     << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
586     cout << "Number of dimensions: " << m_numDim << endl;
587 caltinay 3670
588 caltinay 3691 // write tags
589     if (m_tagMap.size() > 0) {
590     cout << "Tags:" << endl;
591     TagMap::const_iterator it;
592     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
593     cout << " " << setw(5) << it->second << " "
594     << it->first << endl;
595     }
596 caltinay 3670 }
597     }
598    
599     int RipleyDomain::getSystemMatrixTypeId(const int solver,
600     const int preconditioner, const int package, const bool symmetry) const
601     {
602 caltinay 3691 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
603     package, symmetry, m_mpiInfo);
604 caltinay 3670 }
605    
606     int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
607     const int package, const bool symmetry) const
608     {
609 caltinay 3691 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
610     package, symmetry, m_mpiInfo);
611 caltinay 3670 }
612    
613 caltinay 3691 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
614     const escript::FunctionSpace& row_functionspace,
615     const int column_blocksize,
616     const escript::FunctionSpace& column_functionspace,
617     const int type) const
618 caltinay 3670 {
619 caltinay 3691 bool reduceRowOrder=false;
620     bool reduceColOrder=false;
621     // is the domain right?
622     const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
623     if (row_domain!=*this)
624 caltinay 3740 throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
625 caltinay 3691 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
626     if (col_domain!=*this)
627 caltinay 3740 throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
628 caltinay 3691 // is the function space type right?
629 caltinay 3750 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
630 caltinay 3691 reduceRowOrder=true;
631 caltinay 3750 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
632 caltinay 3740 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
633 caltinay 3750 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
634 caltinay 3691 reduceColOrder=true;
635 caltinay 3750 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
636 caltinay 3740 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
637 caltinay 3691
638     // generate matrix
639     Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
640     Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
641     row_blocksize, column_blocksize, FALSE);
642     paso::checkPasoError();
643     Paso_SystemMatrixPattern_free(pattern);
644     escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
645     row_functionspace, column_blocksize, column_functionspace));
646     return sma;
647 caltinay 3670 }
648    
649 caltinay 3748 void RipleyDomain::addPDEToSystem(
650     escript::AbstractSystemMatrix& mat, escript::Data& rhs,
651     const escript::Data& A, const escript::Data& B, const escript::Data& C,
652     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
653     const escript::Data& d, const escript::Data& y,
654     const escript::Data& d_contact, const escript::Data& y_contact,
655     const escript::Data& d_dirac,const escript::Data& y_dirac) const
656     {
657     if (!d_contact.isEmpty() || !y_contact.isEmpty())
658     throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
659    
660     paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
661     if (!sma)
662     throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
663    
664     if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
665     throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
666    
667     //TODO: more input param checks (shape, function space etc)
668    
669     Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
670    
671     if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
672     throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
673    
674     const int numEq=S->logical_row_block_size;
675     const int numComp=S->logical_col_block_size;
676    
677     if (numEq != numComp)
678     throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
679     //TODO: more system matrix checks
680    
681     if (numEq==1)
682     assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
683     else
684     assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
685     }
686    
687 caltinay 3697 void RipleyDomain::setNewX(const escript::Data& arg)
688     {
689     throw RipleyException("setNewX(): Operation not supported");
690     }
691    
692 caltinay 3722 //protected
693 caltinay 3744 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
694 caltinay 3722 {
695     const dim_t numComp = in.getDataPointSize();
696 caltinay 3744 const dim_t dpp = in.getNumDataPointsPerSample();
697 caltinay 3740 out.requireWrite();
698 caltinay 3722 #pragma omp parallel for
699     for (index_t i=0; i<in.getNumSamples(); i++) {
700     const double* src = in.getSampleDataRO(i);
701 caltinay 3744 double* dest = out.getSampleDataRW(i);
702     for (index_t c=0; c<numComp; c++) {
703     double res=0.;
704     for (index_t q=0; q<dpp; q++)
705     res+=src[c+q*numComp];
706     *dest++ = res/dpp;
707     }
708     }
709     }
710    
711     //protected
712     void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
713     {
714     const dim_t numComp = in.getDataPointSize();
715     out.requireWrite();
716     #pragma omp parallel for
717     for (index_t i=0; i<in.getNumSamples(); i++) {
718     const double* src = in.getSampleDataRO(i);
719 caltinay 3722 copy(src, src+numComp, out.getSampleDataRW(i));
720     }
721     }
722    
723     //protected
724     void RipleyDomain::updateTagsInUse(int fsType) const
725     {
726     IndexVector* tagsInUse=NULL;
727     const IndexVector* tags=NULL;
728     switch(fsType) {
729     case Nodes:
730     tags=&m_nodeTags;
731     tagsInUse=&m_nodeTagsInUse;
732     break;
733     case Elements:
734     case ReducedElements:
735     tags=&m_elementTags;
736     tagsInUse=&m_elementTagsInUse;
737     break;
738     case FaceElements:
739     case ReducedFaceElements:
740     tags=&m_faceTags;
741     tagsInUse=&m_faceTagsInUse;
742     break;
743     default:
744     return;
745     }
746    
747     // gather global unique tag values from tags into tagsInUse
748     tagsInUse->clear();
749     index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
750    
751     while (true) {
752     // find smallest value bigger than lastFoundValue
753     minFoundValue = INDEX_T_MAX;
754     #pragma omp parallel private(local_minFoundValue)
755     {
756     local_minFoundValue = minFoundValue;
757     #pragma omp for schedule(static) nowait
758     for (size_t i = 0; i < tags->size(); i++) {
759     const index_t v = (*tags)[i];
760     if ((v > lastFoundValue) && (v < local_minFoundValue))
761     local_minFoundValue = v;
762     }
763     #pragma omp critical
764     {
765     if (local_minFoundValue < minFoundValue)
766     minFoundValue = local_minFoundValue;
767     }
768     }
769     #ifdef ESYS_MPI
770     local_minFoundValue = minFoundValue;
771     MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
772     #endif
773 caltinay 3740
774 caltinay 3722 // if we found a new value add it to the tagsInUse vector
775     if (minFoundValue < INDEX_T_MAX) {
776     tagsInUse->push_back(minFoundValue);
777     lastFoundValue = minFoundValue;
778     } else
779     break;
780     }
781     }
782    
783 caltinay 3756 //protected
784     Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
785     const IndexVector& index, const dim_t M, const dim_t N) const
786     {
787     // paso will manage the memory
788     index_t* indexC = MEMALLOC(index.size(), index_t);
789     index_t* ptrC = MEMALLOC(ptr.size(), index_t);
790     copy(index.begin(), index.end(), indexC);
791     copy(ptr.begin(), ptr.end(), ptrC);
792     return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
793     }
794    
795     //protected
796     Paso_Pattern* RipleyDomain::createMainPattern() const
797     {
798     IndexVector ptr(1,0);
799     IndexVector index;
800    
801     for (index_t i=0; i<getNumDOF(); i++) {
802     // add the DOF itself
803     index.push_back(i);
804     const dim_t num=insertNeighbourNodes(index, i);
805     ptr.push_back(ptr.back()+num+1);
806     }
807    
808     return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
809     }
810    
811     //protected
812     void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
813     const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
814     {
815     IndexVector ptr(1,0);
816     IndexVector index;
817     for (index_t i=0; i<getNumDOF(); i++) {
818     index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
819     ptr.push_back(ptr.back()+colIndices[i].size());
820     }
821    
822     const dim_t M=ptr.size()-1;
823     *colPattern=createPasoPattern(ptr, index, M, N);
824    
825     IndexVector rowPtr(1,0);
826     IndexVector rowIndex;
827     for (dim_t id=0; id<N; id++) {
828     dim_t n=0;
829     for (dim_t i=0; i<M; i++) {
830     for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
831     if (index[j]==id) {
832     rowIndex.push_back(i);
833     n++;
834     break;
835     }
836     }
837     }
838     rowPtr.push_back(rowPtr.back()+n);
839     }
840    
841     // M and N are now reversed
842     *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
843     }
844    
845     //protected
846     void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
847     const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
848     dim_t num_Sol, const vector<double>& array) const
849     {
850     const dim_t numMyCols = mat->pattern->mainPattern->numInput;
851     const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
852    
853     const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
854     const index_t* mainBlock_index = mat->mainBlock->pattern->index;
855     double* mainBlock_val = mat->mainBlock->val;
856     const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
857     const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
858     double* col_coupleBlock_val = mat->col_coupleBlock->val;
859     const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
860     const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
861     double* row_coupleBlock_val = mat->row_coupleBlock->val;
862    
863     for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
864     // down columns of array
865     const dim_t j_Eq = nodes_Eq[k_Eq];
866     const dim_t i_row = j_Eq;
867     //printf("row:%d\n", i_row);
868     // only look at the matrix rows stored on this processor
869     if (i_row < numMyRows) {
870     for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
871     const dim_t i_col = nodes_Sol[k_Sol];
872     if (i_col < numMyCols) {
873     for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
874     if (mainBlock_index[k] == i_col) {
875     mainBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
876     break;
877     }
878     }
879     } else {
880     for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
881     //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;
882     if (col_coupleBlock_index[k] == i_col - numMyCols) {
883     col_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
884     break;
885     }
886     }
887     }
888     }
889     } else {
890     for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
891     // across rows of array
892     const dim_t i_col = nodes_Sol[k_Sol];
893     if (i_col < numMyCols) {
894     for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
895     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
896     {
897     //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;
898     if (row_coupleBlock_index[k] == i_col) {
899     row_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
900     break;
901     }
902     }
903     }
904     }
905     }
906     }
907     }
908    
909 caltinay 3691 //
910     // the methods that follow have to be implemented by the subclasses
911     //
912    
913     string RipleyDomain::getDescription() const
914 caltinay 3670 {
915 caltinay 3691 throw RipleyException("getDescription() not implemented");
916 caltinay 3670 }
917    
918 caltinay 3691 void RipleyDomain::write(const string& filename) const
919 caltinay 3670 {
920 caltinay 3691 throw RipleyException("write() not implemented");
921 caltinay 3670 }
922    
923 caltinay 3691 void RipleyDomain::dump(const string& filename) const
924 caltinay 3670 {
925 caltinay 3691 throw RipleyException("dump() not implemented");
926 caltinay 3670 }
927    
928 caltinay 3691 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
929     {
930     throw RipleyException("borrowSampleReferenceIDs() not implemented");
931     }
932 caltinay 3670
933 caltinay 3691 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
934 caltinay 3670 {
935 caltinay 3691 throw RipleyException("interpolateACross() not implemented");
936 caltinay 3670 }
937    
938 caltinay 3691 bool RipleyDomain::probeInterpolationACross(int fsType_source,
939     const escript::AbstractDomain&, int fsType_target) const
940 caltinay 3670 {
941 caltinay 3691 throw RipleyException("probeInterpolationACross() not implemented");
942 caltinay 3670 }
943    
944 caltinay 3691 void RipleyDomain::setToNormal(escript::Data& normal) const
945 caltinay 3670 {
946 caltinay 3691 throw RipleyException("setToNormal() not implemented");
947 caltinay 3670 }
948    
949 caltinay 3691 void RipleyDomain::setToSize(escript::Data& size) const
950 caltinay 3670 {
951 caltinay 3691 throw RipleyException("setToSize() not implemented");
952 caltinay 3670 }
953    
954 caltinay 3691 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
955 caltinay 3670 {
956 caltinay 3691 throw RipleyException("setToGradient() not implemented");
957 caltinay 3670 }
958    
959 caltinay 3691 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
960 caltinay 3670 {
961 caltinay 3691 throw RipleyException("setToIntegrals() not implemented");
962 caltinay 3670 }
963    
964 caltinay 3691 bool RipleyDomain::ownSample(int fsType, index_t id) const
965 caltinay 3670 {
966 caltinay 3691 throw RipleyException("ownSample() not implemented");
967 caltinay 3670 }
968    
969 caltinay 3691 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
970     const escript::Data& D, const escript::Data& d,
971     const escript::Data& d_dirac, const bool useHRZ) const
972     {
973     throw RipleyException("addPDEToLumpedSystem() not implemented");
974 caltinay 3670 }
975    
976 caltinay 3691 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
977     const escript::Data& Y, const escript::Data& y,
978     const escript::Data& y_contact, const escript::Data& y_dirac) const
979 caltinay 3670 {
980 caltinay 3691 throw RipleyException("addPDEToRHS() not implemented");
981 caltinay 3670 }
982    
983 caltinay 3691 void RipleyDomain::addPDEToTransportProblem(
984     escript::AbstractTransportProblem& tp,
985     escript::Data& source, const escript::Data& M,
986     const escript::Data& A, const escript::Data& B, const escript::Data& C,
987     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
988     const escript::Data& d, const escript::Data& y,
989     const escript::Data& d_contact, const escript::Data& y_contact,
990     const escript::Data& d_dirac, const escript::Data& y_dirac) const
991 caltinay 3670 {
992 caltinay 3691 throw RipleyException("addPDEToTransportProblem() not implemented");
993 caltinay 3670 }
994    
995 caltinay 3691 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
996     const int blocksize, const escript::FunctionSpace& functionspace,
997     const int type) const
998 caltinay 3670 {
999 caltinay 3691 throw RipleyException("newTransportProblem() not implemented");
1000 caltinay 3670 }
1001    
1002 caltinay 3691 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1003     bool reducedColOrder) const
1004 caltinay 3670 {
1005 caltinay 3691 throw RipleyException("getPattern() not implemented");
1006 caltinay 3670 }
1007    
1008 caltinay 3691 dim_t RipleyDomain::getNumDataPointsGlobal() const
1009 caltinay 3670 {
1010 caltinay 3691 throw RipleyException("getNumDataPointsGlobal() not implemented");
1011 caltinay 3670 }
1012    
1013 caltinay 3697 IndexVector RipleyDomain::getNumNodesPerDim() const
1014     {
1015     throw RipleyException("getNumNodesPerDim() not implemented");
1016     }
1017    
1018     IndexVector RipleyDomain::getNumElementsPerDim() const
1019     {
1020     throw RipleyException("getNumElementsPerDim() not implemented");
1021     }
1022    
1023     IndexVector RipleyDomain::getNumFacesPerBoundary() const
1024     {
1025     throw RipleyException("getNumFacesPerBoundary() not implemented");
1026     }
1027    
1028     IndexVector RipleyDomain::getNodeDistribution() const
1029     {
1030     throw RipleyException("getNodeDistribution() not implemented");
1031     }
1032    
1033     pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1034     {
1035     throw RipleyException("getFirstCoordAndSpacing() not implemented");
1036     }
1037    
1038 caltinay 3691 dim_t RipleyDomain::getNumFaceElements() const
1039 caltinay 3670 {
1040 caltinay 3691 throw RipleyException("getNumFaceElements() not implemented");
1041 caltinay 3670 }
1042    
1043 caltinay 3691 dim_t RipleyDomain::getNumElements() const
1044 caltinay 3670 {
1045 caltinay 3691 throw RipleyException("getNumElements() not implemented");
1046 caltinay 3670 }
1047    
1048 caltinay 3691 dim_t RipleyDomain::getNumNodes() const
1049 caltinay 3670 {
1050 caltinay 3691 throw RipleyException("getNumNodes() not implemented");
1051 caltinay 3670 }
1052    
1053 caltinay 3750 dim_t RipleyDomain::getNumDOF() const
1054     {
1055     throw RipleyException("getNumDOF() not implemented");
1056     }
1057    
1058 caltinay 3756 dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1059     {
1060     throw RipleyException("insertNeighbourNodes() not implemented");
1061     }
1062    
1063 caltinay 3691 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1064 caltinay 3670 {
1065 caltinay 3691 throw RipleyException("assembleCoordinates() not implemented");
1066 caltinay 3670 }
1067    
1068 caltinay 3748 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1069     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1070     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1071     const escript::Data& d, const escript::Data& y) const
1072     {
1073     throw RipleyException("assemblePDESingle() not implemented");
1074     }
1075    
1076     void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1077     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1078     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1079     const escript::Data& d, const escript::Data& y) const
1080     {
1081     throw RipleyException("assemblePDESystem() not implemented");
1082     }
1083    
1084 caltinay 3711 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1085 caltinay 3701 {
1086     throw RipleyException("interpolateNodesOnElements() not implemented");
1087     }
1088 caltinay 3691
1089 caltinay 3711 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1090 caltinay 3701 {
1091     throw RipleyException("interpolateNodesOnFaces() not implemented");
1092     }
1093    
1094 caltinay 3750 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1095     {
1096     throw RipleyException("nodesToDOF() not implemented");
1097     }
1098 caltinay 3701
1099 caltinay 3754 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1100     {
1101     throw RipleyException("dofToNodes() not implemented");
1102     }
1103    
1104 caltinay 3670 } // end of namespace ripley
1105    

  ViewVC Help
Powered by ViewVC 1.1.26