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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3697 - (show 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
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 #include <escript/FunctionSpaceFactory.h>
17 #include <pasowrap/SystemMatrixAdapter.h>
18 #include <pasowrap/TransportProblemAdapter.h>
19
20 #include <iomanip>
21
22 using namespace std;
23 using paso::SystemMatrixAdapter;
24 using paso::TransportProblemAdapter;
25
26 namespace ripley {
27
28 escript::Domain_ptr RipleyDomain::loadMesh(const string& filename)
29 {
30 throw RipleyException("loadMesh() not implemented");
31 }
32
33 escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34 {
35 throw RipleyException("readMesh() not implemented");
36 }
37
38 RipleyDomain::RipleyDomain(dim_t dim) :
39 m_numDim(dim),
40 m_status(0)
41 {
42 m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
43 }
44
45 RipleyDomain::~RipleyDomain()
46 {
47 Esys_MPIInfo_free(m_mpiInfo);
48 }
49
50 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
51 {
52 switch (fsType) {
53 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 return true;
63 default:
64 break;
65 }
66 return false;
67 }
68
69 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
70 {
71 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 break;
83 }
84 return "Invalid function space type code";
85 }
86
87 pair<int,int> RipleyDomain::getDataShape(int fsType) const
88 {
89 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
90 switch (fsType) {
91 case Nodes:
92 case DegreesOfFreedom:
93 return pair<int,int>(1, getNumNodes());
94 case Elements:
95 return pair<int,int>(ptsPerSample, getNumElements());
96 case FaceElements:
97 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
98 /*
99 case Points:
100 return pair<int,int>(1, getNumNodes());
101 case ReducedNodes:
102 case ReducedDegreesOfFreedom:
103 return pair<int,int>(1, getNumReducedNodes());
104 case ReducedElements:
105 return pair<int,int>(1, getNumReducedElements());
106 case ReducedFaceElements:
107 return pair<int,int>(1, getNumReducedFaceElements());
108 */
109 default:
110 break;
111 }
112
113 stringstream msg;
114 msg << "Invalid function space type " << fsType << " for "
115 << getDescription();
116 throw RipleyException(msg.str());
117 }
118
119 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
120 {
121 throw RipleyException("getTagFromSampleNo() not implemented");
122 }
123
124 string RipleyDomain::showTagNames() const
125 {
126 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 }
132 return ret.str();
133 }
134
135 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
136 {
137 /*
138 The idea is to use equivalence classes (i.e. types which can be
139 interpolated back and forth):
140 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 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 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 if (hasline[0]==1)
211 resultcode=Points;
212 else if (hasline[1]==1) {
213 if (hasclass[5]==1)
214 resultcode=ReducedElements;
215 else
216 resultcode=Elements;
217 } else if (hasline[2]==1) {
218 if (hasclass[7]==1)
219 resultcode=ReducedFaceElements;
220 else
221 resultcode=FaceElements;
222 } else
223 throw RipleyException("Internal Ripley Error!");
224 } else { // totlines==0
225 if (hasclass[2]==1)
226 // something from class 2
227 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228 else // something from class 1
229 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230 }
231 return true;
232 }
233
234 escript::Data RipleyDomain::getX() const
235 {
236 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 }
266 }
267
268 bool RipleyDomain::isCellOriented(int fsType) const
269 {
270 switch(fsType) {
271 case Nodes:
272 case DegreesOfFreedom:
273 case ReducedDegreesOfFreedom:
274 return false;
275 case Elements:
276 case FaceElements:
277 case Points:
278 case ReducedElements:
279 case ReducedFaceElements:
280 return true;
281 default:
282 break;
283 }
284 stringstream msg;
285 msg << "Illegal function space type " << fsType << " on " << getDescription();
286 throw RipleyException(msg.str());
287 }
288
289 void RipleyDomain::Print_Mesh_Info(const bool full) const
290 {
291 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
295 // 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 }
304 }
305
306 int RipleyDomain::getSystemMatrixTypeId(const int solver,
307 const int preconditioner, const int package, const bool symmetry) const
308 {
309 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
310 package, symmetry, m_mpiInfo);
311 }
312
313 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
314 const int package, const bool symmetry) const
315 {
316 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
317 package, symmetry, m_mpiInfo);
318 }
319
320 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 {
326 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 }
355
356 void RipleyDomain::setNewX(const escript::Data& arg)
357 {
358 throw RipleyException("setNewX(): Operation not supported");
359 }
360
361 //
362 // the methods that follow have to be implemented by the subclasses
363 //
364
365 string RipleyDomain::getDescription() const
366 {
367 throw RipleyException("getDescription() not implemented");
368 }
369
370 bool RipleyDomain::operator==(const AbstractDomain& other) const
371 {
372 throw RipleyException("operator==() not implemented");
373 }
374
375 void RipleyDomain::write(const string& filename) const
376 {
377 throw RipleyException("write() not implemented");
378 }
379
380 void RipleyDomain::dump(const string& filename) const
381 {
382 throw RipleyException("dump() not implemented");
383 }
384
385 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
386 {
387 throw RipleyException("borrowSampleReferenceIDs() not implemented");
388 }
389
390 void RipleyDomain::interpolateOnDomain(escript::Data& target,
391 const escript::Data& in) const
392 {
393 throw RipleyException("interpolateOnDomain() not implemented");
394 }
395
396 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
397 int fsType_target) const
398 {
399 throw RipleyException("probeInterpolationOnDomain() not implemented");
400 }
401
402 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
403 {
404 throw RipleyException("interpolateACross() not implemented");
405 }
406
407 bool RipleyDomain::probeInterpolationACross(int fsType_source,
408 const escript::AbstractDomain&, int fsType_target) const
409 {
410 throw RipleyException("probeInterpolationACross() not implemented");
411 }
412
413 void RipleyDomain::setToNormal(escript::Data& normal) const
414 {
415 throw RipleyException("setToNormal() not implemented");
416 }
417
418 void RipleyDomain::setToSize(escript::Data& size) const
419 {
420 throw RipleyException("setToSize() not implemented");
421 }
422
423 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
424 {
425 throw RipleyException("setToGradient() not implemented");
426 }
427
428 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
429 {
430 throw RipleyException("setToIntegrals() not implemented");
431 }
432
433 bool RipleyDomain::ownSample(int fsType, index_t id) const
434 {
435 throw RipleyException("ownSample() not implemented");
436 }
437
438 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const
439 {
440 throw RipleyException("setTags() not implemented");
441 }
442
443 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
444 {
445 throw RipleyException("getNumberOfTagsInUse() not implemented");
446 }
447
448 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
449 {
450 throw RipleyException("borrowListOfTagsInUse() not implemented");
451 }
452
453 bool RipleyDomain::canTag(int fsType) const
454 {
455 throw RipleyException("canTag() not implemented");
456 }
457
458
459 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
470 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 }
476
477 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 {
481 throw RipleyException("addPDEToRHS() not implemented");
482 }
483
484 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 {
493 throw RipleyException("addPDEToTransportProblem() not implemented");
494 }
495
496 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
497 const int blocksize, const escript::FunctionSpace& functionspace,
498 const int type) const
499 {
500 throw RipleyException("newTransportProblem() not implemented");
501 }
502
503 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
504 bool reducedColOrder) const
505 {
506 throw RipleyException("getPattern() not implemented");
507 }
508
509 dim_t RipleyDomain::getNumDataPointsGlobal() const
510 {
511 throw RipleyException("getNumDataPointsGlobal() not implemented");
512 }
513
514 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 dim_t RipleyDomain::getNumFaceElements() const
540 {
541 throw RipleyException("getNumFaceElements() not implemented");
542 }
543
544 dim_t RipleyDomain::getNumElements() const
545 {
546 throw RipleyException("getNumElements() not implemented");
547 }
548
549 dim_t RipleyDomain::getNumNodes() const
550 {
551 throw RipleyException("getNumNodes() not implemented");
552 }
553
554 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
555 {
556 throw RipleyException("assembleCoordinates() not implemented");
557 }
558
559
560 } // end of namespace ripley
561

  ViewVC Help
Powered by ViewVC 1.1.26