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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3691 - (hide annotations)
Wed Nov 23 23:07:37 2011 UTC (7 years, 9 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 17106 byte(s)
Next iteration

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

  ViewVC Help
Powered by ViewVC 1.1.26