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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3750 - (hide annotations)
Fri Dec 23 01:20:34 2011 UTC (7 years, 9 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 30405 byte(s)
Checkpoint - reinstated DOF FS and tried to fix couple blocks...

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

  ViewVC Help
Powered by ViewVC 1.1.26