/[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 3697 - (hide annotations)
Mon Nov 28 04:52:00 2011 UTC (7 years, 10 months ago) by caltinay
File size: 17659 byte(s)
[x] Rectangle node id's and weipa compatibility
[ ] Brick...

1 caltinay 3670
2     /*******************************************************
3     *
4     * Copyright (c) 2003-2011 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14     #include <ripley/RipleyDomain.h>
15     #include <escript/DataFactory.h>
16 caltinay 3691 #include <escript/FunctionSpaceFactory.h>
17     #include <pasowrap/SystemMatrixAdapter.h>
18     #include <pasowrap/TransportProblemAdapter.h>
19 caltinay 3670
20     #include <iomanip>
21    
22     using namespace std;
23 caltinay 3691 using paso::SystemMatrixAdapter;
24     using paso::TransportProblemAdapter;
25 caltinay 3670
26     namespace ripley {
27    
28 caltinay 3691 escript::Domain_ptr RipleyDomain::loadMesh(const string& filename)
29 caltinay 3670 {
30 caltinay 3691 throw RipleyException("loadMesh() not implemented");
31 caltinay 3670 }
32    
33 caltinay 3691 escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34 caltinay 3670 {
35 caltinay 3691 throw RipleyException("readMesh() not implemented");
36 caltinay 3670 }
37    
38 caltinay 3691 RipleyDomain::RipleyDomain(dim_t dim) :
39     m_numDim(dim),
40     m_status(0)
41 caltinay 3670 {
42 caltinay 3691 m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
43 caltinay 3670 }
44    
45 caltinay 3691 RipleyDomain::~RipleyDomain()
46 caltinay 3670 {
47 caltinay 3691 Esys_MPIInfo_free(m_mpiInfo);
48 caltinay 3670 }
49    
50 caltinay 3691 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
51 caltinay 3670 {
52 caltinay 3691 switch (fsType) {
53 caltinay 3670 case DegreesOfFreedom:
54     case ReducedDegreesOfFreedom:
55     case Nodes:
56     case ReducedNodes:
57     case Elements:
58     case ReducedElements:
59     case FaceElements:
60     case ReducedFaceElements:
61     case Points:
62 caltinay 3691 return true;
63     default:
64 caltinay 3670 break;
65     }
66 caltinay 3691 return false;
67 caltinay 3670 }
68    
69 caltinay 3691 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
70 caltinay 3670 {
71 caltinay 3691 switch (fsType) {
72     case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
73     case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
74     case Nodes: return "Ripley_Nodes";
75     case ReducedNodes: return "Ripley_Reduced_Nodes";
76     case Elements: return "Ripley_Elements";
77     case ReducedElements: return "Ripley_Reduced_Elements";
78     case FaceElements: return "Ripley_Face_Elements";
79     case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";
80     case Points: return "Ripley_Points";
81     default:
82 caltinay 3670 break;
83     }
84 caltinay 3691 return "Invalid function space type code";
85 caltinay 3670 }
86    
87 caltinay 3691 pair<int,int> RipleyDomain::getDataShape(int fsType) const
88 caltinay 3670 {
89 caltinay 3691 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
90     switch (fsType) {
91 caltinay 3670 case Nodes:
92     case DegreesOfFreedom:
93 caltinay 3691 return pair<int,int>(1, getNumNodes());
94 caltinay 3670 case Elements:
95 caltinay 3691 return pair<int,int>(ptsPerSample, getNumElements());
96 caltinay 3670 case FaceElements:
97 caltinay 3691 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
98     /*
99 caltinay 3697 case Points:
100     return pair<int,int>(1, getNumNodes());
101 caltinay 3670 case ReducedNodes:
102 caltinay 3697 case ReducedDegreesOfFreedom:
103 caltinay 3691 return pair<int,int>(1, getNumReducedNodes());
104 caltinay 3670 case ReducedElements:
105 caltinay 3691 return pair<int,int>(1, getNumReducedElements());
106 caltinay 3670 case ReducedFaceElements:
107 caltinay 3691 return pair<int,int>(1, getNumReducedFaceElements());
108     */
109 caltinay 3670 default:
110     break;
111     }
112    
113 caltinay 3691 stringstream msg;
114     msg << "Invalid function space type " << fsType << " for "
115     << getDescription();
116     throw RipleyException(msg.str());
117 caltinay 3670 }
118    
119 caltinay 3691 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
120 caltinay 3670 {
121 caltinay 3691 throw RipleyException("getTagFromSampleNo() not implemented");
122 caltinay 3670 }
123    
124 caltinay 3691 string RipleyDomain::showTagNames() const
125 caltinay 3670 {
126 caltinay 3691 stringstream ret;
127     TagMap::const_iterator it;
128     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
129     if (it!=m_tagMap.begin()) ret << ", ";
130     ret << it->first;
131 caltinay 3670 }
132 caltinay 3691 return ret.str();
133 caltinay 3670 }
134    
135     bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
136     {
137 caltinay 3691 /*
138     The idea is to use equivalence classes (i.e. types which can be
139     interpolated back and forth):
140 caltinay 3670 class 1: DOF <-> Nodes
141     class 2: ReducedDOF <-> ReducedNodes
142     class 3: Points
143     class 4: Elements
144     class 5: ReducedElements
145     class 6: FaceElements
146     class 7: ReducedFaceElements
147     class 8: ContactElementZero <-> ContactElementOne
148     class 9: ReducedContactElementZero <-> ReducedContactElementOne
149    
150 caltinay 3691 There is also a set of lines. Interpolation is possible down a line but not
151     between lines.
152     class 1 and 2 belong to all lines so aren't considered.
153 caltinay 3670 line 0: class 3
154     line 1: class 4,5
155     line 2: class 6,7
156     line 3: class 8,9
157    
158     For classes with multiple members (eg class 2) we have vars to record if
159     there is at least one instance. e.g. hasnodes is true if we have at least
160     one instance of Nodes.
161     */
162     if (fs.empty())
163     return false;
164     vector<int> hasclass(10);
165     vector<int> hasline(4);
166     bool hasnodes=false;
167     bool hasrednodes=false;
168     for (int i=0;i<fs.size();++i) {
169     switch (fs[i]) {
170     case Nodes: hasnodes=true; // no break is deliberate
171     case DegreesOfFreedom:
172     hasclass[1]=1;
173     break;
174     case ReducedNodes: hasrednodes=true; // no break is deliberate
175     case ReducedDegreesOfFreedom:
176     hasclass[2]=1;
177     break;
178     case Points:
179     hasline[0]=1;
180     hasclass[3]=1;
181     break;
182     case Elements:
183     hasclass[4]=1;
184     hasline[1]=1;
185     break;
186     case ReducedElements:
187     hasclass[5]=1;
188     hasline[1]=1;
189     break;
190     case FaceElements:
191     hasclass[6]=1;
192     hasline[2]=1;
193     break;
194     case ReducedFaceElements:
195     hasclass[7]=1;
196     hasline[2]=1;
197     break;
198     default:
199     return false;
200     }
201     }
202     int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
203    
204     // fail if we have more than one leaf group
205     // = there are at least two branches we can't interpolate between
206     if (totlines>1)
207     return false;
208     else if (totlines==1) {
209     // we have points
210 caltinay 3691 if (hasline[0]==1)
211 caltinay 3670 resultcode=Points;
212 caltinay 3691 else if (hasline[1]==1) {
213     if (hasclass[5]==1)
214 caltinay 3670 resultcode=ReducedElements;
215 caltinay 3691 else
216 caltinay 3670 resultcode=Elements;
217     } else if (hasline[2]==1) {
218 caltinay 3691 if (hasclass[7]==1)
219 caltinay 3670 resultcode=ReducedFaceElements;
220 caltinay 3691 else
221 caltinay 3670 resultcode=FaceElements;
222 caltinay 3691 } else
223     throw RipleyException("Internal Ripley Error!");
224 caltinay 3670 } else { // totlines==0
225 caltinay 3691 if (hasclass[2]==1)
226 caltinay 3670 // something from class 2
227     resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228 caltinay 3691 else // something from class 1
229 caltinay 3670 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230     }
231     return true;
232     }
233    
234 caltinay 3691 escript::Data RipleyDomain::getX() const
235 caltinay 3670 {
236 caltinay 3691 return escript::continuousFunction(*this).getX();
237     }
238    
239     escript::Data RipleyDomain::getNormal() const
240     {
241     return escript::functionOnBoundary(*this).getNormal();
242     }
243    
244     escript::Data RipleyDomain::getSize() const
245     {
246     return escript::function(*this).getSize();
247     }
248    
249     void RipleyDomain::setToX(escript::Data& arg) const
250     {
251     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
252     *(arg.getFunctionSpace().getDomain()));
253     if (argDomain != *this)
254     throw RipleyException("setToX: Illegal domain of data point locations");
255     if (!arg.isExpanded())
256     throw RipleyException("setToX: Expanded Data object expected");
257    
258     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
259     assembleCoordinates(arg);
260     } else {
261     // interpolate the result
262     escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
263     assembleCoordinates(contData);
264     interpolateOnDomain(arg, contData);
265 caltinay 3670 }
266 caltinay 3691 }
267 caltinay 3670
268 caltinay 3691 bool RipleyDomain::isCellOriented(int fsType) const
269     {
270     switch(fsType) {
271 caltinay 3670 case Nodes:
272     case DegreesOfFreedom:
273     case ReducedDegreesOfFreedom:
274 caltinay 3691 return false;
275 caltinay 3670 case Elements:
276 caltinay 3691 case FaceElements:
277     case Points:
278 caltinay 3670 case ReducedElements:
279     case ReducedFaceElements:
280 caltinay 3691 return true;
281 caltinay 3670 default:
282 caltinay 3691 break;
283 caltinay 3670 }
284 caltinay 3691 stringstream msg;
285     msg << "Illegal function space type " << fsType << " on " << getDescription();
286     throw RipleyException(msg.str());
287 caltinay 3670 }
288    
289 caltinay 3691 void RipleyDomain::Print_Mesh_Info(const bool full) const
290 caltinay 3670 {
291 caltinay 3691 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
292     << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
293     cout << "Number of dimensions: " << m_numDim << endl;
294 caltinay 3670
295 caltinay 3691 // write tags
296     if (m_tagMap.size() > 0) {
297     cout << "Tags:" << endl;
298     TagMap::const_iterator it;
299     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
300     cout << " " << setw(5) << it->second << " "
301     << it->first << endl;
302     }
303 caltinay 3670 }
304     }
305    
306     int RipleyDomain::getSystemMatrixTypeId(const int solver,
307     const int preconditioner, const int package, const bool symmetry) const
308     {
309 caltinay 3691 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
310     package, symmetry, m_mpiInfo);
311 caltinay 3670 }
312    
313     int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
314     const int package, const bool symmetry) const
315     {
316 caltinay 3691 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
317     package, symmetry, m_mpiInfo);
318 caltinay 3670 }
319    
320 caltinay 3691 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
321     const escript::FunctionSpace& row_functionspace,
322     const int column_blocksize,
323     const escript::FunctionSpace& column_functionspace,
324     const int type) const
325 caltinay 3670 {
326 caltinay 3691 bool reduceRowOrder=false;
327     bool reduceColOrder=false;
328     // is the domain right?
329     const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
330     if (row_domain!=*this)
331     throw RipleyException("Domain of row function space does not match the domain of matrix generator");
332     const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
333     if (col_domain!=*this)
334     throw RipleyException("Domain of column function space does not match the domain of matrix generator");
335     // is the function space type right?
336     if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
337     reduceRowOrder=true;
338     else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
339     throw RipleyException("Illegal function space type for system matrix rows");
340     if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
341     reduceColOrder=true;
342     else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
343     throw RipleyException("Illegal function space type for system matrix columns");
344    
345     // generate matrix
346     Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
347     Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
348     row_blocksize, column_blocksize, FALSE);
349     paso::checkPasoError();
350     Paso_SystemMatrixPattern_free(pattern);
351     escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
352     row_functionspace, column_blocksize, column_functionspace));
353     return sma;
354 caltinay 3670 }
355    
356 caltinay 3697 void RipleyDomain::setNewX(const escript::Data& arg)
357     {
358     throw RipleyException("setNewX(): Operation not supported");
359     }
360    
361 caltinay 3691 //
362     // the methods that follow have to be implemented by the subclasses
363     //
364    
365     string RipleyDomain::getDescription() const
366 caltinay 3670 {
367 caltinay 3691 throw RipleyException("getDescription() not implemented");
368 caltinay 3670 }
369    
370 caltinay 3691 bool RipleyDomain::operator==(const AbstractDomain& other) const
371 caltinay 3670 {
372 caltinay 3691 throw RipleyException("operator==() not implemented");
373 caltinay 3670 }
374    
375 caltinay 3691 void RipleyDomain::write(const string& filename) const
376 caltinay 3670 {
377 caltinay 3691 throw RipleyException("write() not implemented");
378 caltinay 3670 }
379    
380 caltinay 3691 void RipleyDomain::dump(const string& filename) const
381 caltinay 3670 {
382 caltinay 3691 throw RipleyException("dump() not implemented");
383 caltinay 3670 }
384    
385 caltinay 3691 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
386     {
387     throw RipleyException("borrowSampleReferenceIDs() not implemented");
388     }
389 caltinay 3670
390 caltinay 3691 void RipleyDomain::interpolateOnDomain(escript::Data& target,
391     const escript::Data& in) const
392 caltinay 3670 {
393 caltinay 3691 throw RipleyException("interpolateOnDomain() not implemented");
394 caltinay 3670 }
395    
396 caltinay 3691 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
397     int fsType_target) const
398 caltinay 3670 {
399 caltinay 3691 throw RipleyException("probeInterpolationOnDomain() not implemented");
400 caltinay 3670 }
401    
402 caltinay 3691 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
403 caltinay 3670 {
404 caltinay 3691 throw RipleyException("interpolateACross() not implemented");
405 caltinay 3670 }
406    
407 caltinay 3691 bool RipleyDomain::probeInterpolationACross(int fsType_source,
408     const escript::AbstractDomain&, int fsType_target) const
409 caltinay 3670 {
410 caltinay 3691 throw RipleyException("probeInterpolationACross() not implemented");
411 caltinay 3670 }
412    
413 caltinay 3691 void RipleyDomain::setToNormal(escript::Data& normal) const
414 caltinay 3670 {
415 caltinay 3691 throw RipleyException("setToNormal() not implemented");
416 caltinay 3670 }
417    
418 caltinay 3691 void RipleyDomain::setToSize(escript::Data& size) const
419 caltinay 3670 {
420 caltinay 3691 throw RipleyException("setToSize() not implemented");
421 caltinay 3670 }
422    
423 caltinay 3691 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
424 caltinay 3670 {
425 caltinay 3691 throw RipleyException("setToGradient() not implemented");
426 caltinay 3670 }
427    
428 caltinay 3691 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
429 caltinay 3670 {
430 caltinay 3691 throw RipleyException("setToIntegrals() not implemented");
431 caltinay 3670 }
432    
433 caltinay 3691 bool RipleyDomain::ownSample(int fsType, index_t id) const
434 caltinay 3670 {
435 caltinay 3691 throw RipleyException("ownSample() not implemented");
436 caltinay 3670 }
437    
438 caltinay 3691 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const
439 caltinay 3670 {
440 caltinay 3691 throw RipleyException("setTags() not implemented");
441 caltinay 3670 }
442    
443 caltinay 3691 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
444 caltinay 3670 {
445 caltinay 3691 throw RipleyException("getNumberOfTagsInUse() not implemented");
446 caltinay 3670 }
447    
448 caltinay 3691 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
449 caltinay 3670 {
450 caltinay 3691 throw RipleyException("borrowListOfTagsInUse() not implemented");
451 caltinay 3670 }
452    
453 caltinay 3691 bool RipleyDomain::canTag(int fsType) const
454 caltinay 3670 {
455 caltinay 3691 throw RipleyException("canTag() not implemented");
456     }
457 caltinay 3670
458    
459 caltinay 3691 void RipleyDomain::addPDEToSystem(
460     escript::AbstractSystemMatrix& mat, escript::Data& rhs,
461     const escript::Data& A, const escript::Data& B, const escript::Data& C,
462     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
463     const escript::Data& d, const escript::Data& y,
464     const escript::Data& d_contact, const escript::Data& y_contact,
465     const escript::Data& d_dirac,const escript::Data& y_dirac) const
466     {
467     throw RipleyException("addPDEToSystem() not implemented");
468     }
469 caltinay 3670
470 caltinay 3691 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
471     const escript::Data& D, const escript::Data& d,
472     const escript::Data& d_dirac, const bool useHRZ) const
473     {
474     throw RipleyException("addPDEToLumpedSystem() not implemented");
475 caltinay 3670 }
476    
477 caltinay 3691 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
478     const escript::Data& Y, const escript::Data& y,
479     const escript::Data& y_contact, const escript::Data& y_dirac) const
480 caltinay 3670 {
481 caltinay 3691 throw RipleyException("addPDEToRHS() not implemented");
482 caltinay 3670 }
483    
484 caltinay 3691 void RipleyDomain::addPDEToTransportProblem(
485     escript::AbstractTransportProblem& tp,
486     escript::Data& source, const escript::Data& M,
487     const escript::Data& A, const escript::Data& B, const escript::Data& C,
488     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
489     const escript::Data& d, const escript::Data& y,
490     const escript::Data& d_contact, const escript::Data& y_contact,
491     const escript::Data& d_dirac, const escript::Data& y_dirac) const
492 caltinay 3670 {
493 caltinay 3691 throw RipleyException("addPDEToTransportProblem() not implemented");
494 caltinay 3670 }
495    
496 caltinay 3691 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
497     const int blocksize, const escript::FunctionSpace& functionspace,
498     const int type) const
499 caltinay 3670 {
500 caltinay 3691 throw RipleyException("newTransportProblem() not implemented");
501 caltinay 3670 }
502    
503 caltinay 3691 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
504     bool reducedColOrder) const
505 caltinay 3670 {
506 caltinay 3691 throw RipleyException("getPattern() not implemented");
507 caltinay 3670 }
508    
509 caltinay 3691 dim_t RipleyDomain::getNumDataPointsGlobal() const
510 caltinay 3670 {
511 caltinay 3691 throw RipleyException("getNumDataPointsGlobal() not implemented");
512 caltinay 3670 }
513    
514 caltinay 3697 IndexVector RipleyDomain::getNumNodesPerDim() const
515     {
516     throw RipleyException("getNumNodesPerDim() not implemented");
517     }
518    
519     IndexVector RipleyDomain::getNumElementsPerDim() const
520     {
521     throw RipleyException("getNumElementsPerDim() not implemented");
522     }
523    
524     IndexVector RipleyDomain::getNumFacesPerBoundary() const
525     {
526     throw RipleyException("getNumFacesPerBoundary() not implemented");
527     }
528    
529     IndexVector RipleyDomain::getNodeDistribution() const
530     {
531     throw RipleyException("getNodeDistribution() not implemented");
532     }
533    
534     pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
535     {
536     throw RipleyException("getFirstCoordAndSpacing() not implemented");
537     }
538    
539 caltinay 3691 dim_t RipleyDomain::getNumFaceElements() const
540 caltinay 3670 {
541 caltinay 3691 throw RipleyException("getNumFaceElements() not implemented");
542 caltinay 3670 }
543    
544 caltinay 3691 dim_t RipleyDomain::getNumElements() const
545 caltinay 3670 {
546 caltinay 3691 throw RipleyException("getNumElements() not implemented");
547 caltinay 3670 }
548    
549 caltinay 3691 dim_t RipleyDomain::getNumNodes() const
550 caltinay 3670 {
551 caltinay 3691 throw RipleyException("getNumNodes() not implemented");
552 caltinay 3670 }
553    
554 caltinay 3691 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
555 caltinay 3670 {
556 caltinay 3691 throw RipleyException("assembleCoordinates() not implemented");
557 caltinay 3670 }
558    
559 caltinay 3691
560 caltinay 3670 } // end of namespace ripley
561    

  ViewVC Help
Powered by ViewVC 1.1.26