/[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 3713 - (show annotations)
Tue Dec 6 04:43:29 2011 UTC (7 years, 10 months ago) by caltinay
File size: 19801 byte(s)
Integrals in 2D and 3D for rectangle & brick elements & face elements in
regular and reduced orders.

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 ReducedNodes:
93 case DegreesOfFreedom:
94 return pair<int,int>(1, getNumNodes());
95 case Elements:
96 return pair<int,int>(ptsPerSample, getNumElements());
97 case FaceElements:
98 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
99 case ReducedElements:
100 return pair<int,int>(1, getNumElements());
101 case ReducedFaceElements:
102 return pair<int,int>(1, getNumFaceElements());
103 /*
104 case Points:
105 return pair<int,int>(1, getNumPoints());
106 case ReducedDegreesOfFreedom:
107 return pair<int,int>(1, getNumReducedNodes());
108 */
109 default:
110 break;
111 }
112
113 stringstream msg;
114 msg << "getDataShape(): Unsupported function space type "
115 << functionSpaceTypeAsString(fsType) << " for " << getDescription();
116 throw RipleyException(msg.str());
117 }
118
119 string RipleyDomain::showTagNames() const
120 {
121 stringstream ret;
122 TagMap::const_iterator it;
123 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
124 if (it!=m_tagMap.begin()) ret << ", ";
125 ret << it->first;
126 }
127 return ret.str();
128 }
129
130 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
131 {
132 /*
133 The idea is to use equivalence classes (i.e. types which can be
134 interpolated back and forth):
135 class 1: DOF <-> Nodes
136 class 2: ReducedDOF <-> ReducedNodes
137 class 3: Points
138 class 4: Elements
139 class 5: ReducedElements
140 class 6: FaceElements
141 class 7: ReducedFaceElements
142 class 8: ContactElementZero <-> ContactElementOne
143 class 9: ReducedContactElementZero <-> ReducedContactElementOne
144
145 There is also a set of lines. Interpolation is possible down a line but not
146 between lines.
147 class 1 and 2 belong to all lines so aren't considered.
148 line 0: class 3
149 line 1: class 4,5
150 line 2: class 6,7
151 line 3: class 8,9
152
153 For classes with multiple members (eg class 2) we have vars to record if
154 there is at least one instance. e.g. hasnodes is true if we have at least
155 one instance of Nodes.
156 */
157 if (fs.empty())
158 return false;
159 vector<int> hasclass(10);
160 vector<int> hasline(4);
161 bool hasnodes=false;
162 bool hasrednodes=false;
163 for (int i=0;i<fs.size();++i) {
164 switch (fs[i]) {
165 case Nodes: hasnodes=true; // no break is deliberate
166 case DegreesOfFreedom:
167 hasclass[1]=1;
168 break;
169 case ReducedNodes: hasrednodes=true; // no break is deliberate
170 case ReducedDegreesOfFreedom:
171 hasclass[2]=1;
172 break;
173 case Points:
174 hasline[0]=1;
175 hasclass[3]=1;
176 break;
177 case Elements:
178 hasclass[4]=1;
179 hasline[1]=1;
180 break;
181 case ReducedElements:
182 hasclass[5]=1;
183 hasline[1]=1;
184 break;
185 case FaceElements:
186 hasclass[6]=1;
187 hasline[2]=1;
188 break;
189 case ReducedFaceElements:
190 hasclass[7]=1;
191 hasline[2]=1;
192 break;
193 default:
194 return false;
195 }
196 }
197 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
198
199 // fail if we have more than one leaf group
200 // = there are at least two branches we can't interpolate between
201 if (totlines>1)
202 return false;
203 else if (totlines==1) {
204 // we have points
205 if (hasline[0]==1)
206 resultcode=Points;
207 else if (hasline[1]==1) {
208 if (hasclass[5]==1)
209 resultcode=ReducedElements;
210 else
211 resultcode=Elements;
212 } else if (hasline[2]==1) {
213 if (hasclass[7]==1)
214 resultcode=ReducedFaceElements;
215 else
216 resultcode=FaceElements;
217 } else
218 throw RipleyException("Internal Ripley Error!");
219 } else { // totlines==0
220 if (hasclass[2]==1)
221 // something from class 2
222 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
223 else // something from class 1
224 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
225 }
226 return true;
227 }
228
229 void RipleyDomain::interpolateOnDomain(escript::Data& target,
230 const escript::Data& in) const
231 {
232 const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
233 const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
234 if (inDomain != *this)
235 throw RipleyException("Illegal domain of interpolant");
236 if (targetDomain != *this)
237 throw RipleyException("Illegal domain of interpolation target");
238
239 stringstream msg;
240 msg << "interpolateOnDomain() not implemented for function space "
241 << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
242 << " -> "
243 << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
244
245 switch (in.getFunctionSpace().getTypeCode()) {
246 case Nodes:
247 case DegreesOfFreedom:
248 switch (target.getFunctionSpace().getTypeCode()) {
249 case DegreesOfFreedom:
250 // FIXME!
251 target=in;
252 break;
253
254 case Elements:
255 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
256 break;
257
258 case ReducedElements:
259 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
260 break;
261
262 case FaceElements:
263 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
264 break;
265
266 case ReducedFaceElements:
267 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
268 break;
269
270 default:
271 throw RipleyException(msg.str());
272 }
273 break;
274 default:
275 throw RipleyException(msg.str());
276 }
277 }
278
279 escript::Data RipleyDomain::getX() const
280 {
281 return escript::continuousFunction(*this).getX();
282 }
283
284 escript::Data RipleyDomain::getNormal() const
285 {
286 return escript::functionOnBoundary(*this).getNormal();
287 }
288
289 escript::Data RipleyDomain::getSize() const
290 {
291 return escript::function(*this).getSize();
292 }
293
294 void RipleyDomain::setToX(escript::Data& arg) const
295 {
296 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
297 *(arg.getFunctionSpace().getDomain()));
298 if (argDomain != *this)
299 throw RipleyException("setToX: Illegal domain of data point locations");
300 if (!arg.isExpanded())
301 throw RipleyException("setToX: Expanded Data object expected");
302
303 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
304 assembleCoordinates(arg);
305 } else {
306 // interpolate the result
307 escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
308 assembleCoordinates(contData);
309 interpolateOnDomain(arg, contData);
310 }
311 }
312
313 bool RipleyDomain::isCellOriented(int fsType) const
314 {
315 switch(fsType) {
316 case Nodes:
317 case DegreesOfFreedom:
318 case ReducedDegreesOfFreedom:
319 return false;
320 case Elements:
321 case FaceElements:
322 case Points:
323 case ReducedElements:
324 case ReducedFaceElements:
325 return true;
326 default:
327 break;
328 }
329 stringstream msg;
330 msg << "Illegal function space type " << fsType << " on " << getDescription();
331 throw RipleyException(msg.str());
332 }
333
334 void RipleyDomain::Print_Mesh_Info(const bool full) const
335 {
336 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
337 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
338 cout << "Number of dimensions: " << m_numDim << endl;
339
340 // write tags
341 if (m_tagMap.size() > 0) {
342 cout << "Tags:" << endl;
343 TagMap::const_iterator it;
344 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
345 cout << " " << setw(5) << it->second << " "
346 << it->first << endl;
347 }
348 }
349 }
350
351 int RipleyDomain::getSystemMatrixTypeId(const int solver,
352 const int preconditioner, const int package, const bool symmetry) const
353 {
354 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
355 package, symmetry, m_mpiInfo);
356 }
357
358 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
359 const int package, const bool symmetry) const
360 {
361 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
362 package, symmetry, m_mpiInfo);
363 }
364
365 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
366 const escript::FunctionSpace& row_functionspace,
367 const int column_blocksize,
368 const escript::FunctionSpace& column_functionspace,
369 const int type) const
370 {
371 bool reduceRowOrder=false;
372 bool reduceColOrder=false;
373 // is the domain right?
374 const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
375 if (row_domain!=*this)
376 throw RipleyException("Domain of row function space does not match the domain of matrix generator");
377 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
378 if (col_domain!=*this)
379 throw RipleyException("Domain of column function space does not match the domain of matrix generator");
380 // is the function space type right?
381 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
382 reduceRowOrder=true;
383 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
384 throw RipleyException("Illegal function space type for system matrix rows");
385 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
386 reduceColOrder=true;
387 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
388 throw RipleyException("Illegal function space type for system matrix columns");
389
390 // generate matrix
391 Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
392 Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
393 row_blocksize, column_blocksize, FALSE);
394 paso::checkPasoError();
395 Paso_SystemMatrixPattern_free(pattern);
396 escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
397 row_functionspace, column_blocksize, column_functionspace));
398 return sma;
399 }
400
401 void RipleyDomain::setNewX(const escript::Data& arg)
402 {
403 throw RipleyException("setNewX(): Operation not supported");
404 }
405
406 //
407 // the methods that follow have to be implemented by the subclasses
408 //
409
410 string RipleyDomain::getDescription() const
411 {
412 throw RipleyException("getDescription() not implemented");
413 }
414
415 bool RipleyDomain::operator==(const AbstractDomain& other) const
416 {
417 throw RipleyException("operator==() not implemented");
418 }
419
420 void RipleyDomain::write(const string& filename) const
421 {
422 throw RipleyException("write() not implemented");
423 }
424
425 void RipleyDomain::dump(const string& filename) const
426 {
427 throw RipleyException("dump() not implemented");
428 }
429
430 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
431 {
432 throw RipleyException("borrowSampleReferenceIDs() not implemented");
433 }
434
435 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
436 int fsType_target) const
437 {
438 throw RipleyException("probeInterpolationOnDomain() not implemented");
439 }
440
441 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
442 {
443 throw RipleyException("interpolateACross() not implemented");
444 }
445
446 bool RipleyDomain::probeInterpolationACross(int fsType_source,
447 const escript::AbstractDomain&, int fsType_target) const
448 {
449 throw RipleyException("probeInterpolationACross() not implemented");
450 }
451
452 void RipleyDomain::setToNormal(escript::Data& normal) const
453 {
454 throw RipleyException("setToNormal() not implemented");
455 }
456
457 void RipleyDomain::setToSize(escript::Data& size) const
458 {
459 throw RipleyException("setToSize() not implemented");
460 }
461
462 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
463 {
464 throw RipleyException("setToGradient() not implemented");
465 }
466
467 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
468 {
469 throw RipleyException("setToIntegrals() not implemented");
470 }
471
472 bool RipleyDomain::ownSample(int fsType, index_t id) const
473 {
474 throw RipleyException("ownSample() not implemented");
475 }
476
477 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
478 {
479 throw RipleyException("getTagFromSampleNo() not implemented");
480 }
481
482 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const
483 {
484 throw RipleyException("setTags() not implemented");
485 }
486
487 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
488 {
489 throw RipleyException("getNumberOfTagsInUse() not implemented");
490 }
491
492 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
493 {
494 throw RipleyException("borrowListOfTagsInUse() not implemented");
495 }
496
497 bool RipleyDomain::canTag(int fsType) const
498 {
499 return false;
500 //throw RipleyException("canTag() not implemented");
501 }
502
503 void RipleyDomain::addPDEToSystem(
504 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
505 const escript::Data& A, const escript::Data& B, const escript::Data& C,
506 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
507 const escript::Data& d, const escript::Data& y,
508 const escript::Data& d_contact, const escript::Data& y_contact,
509 const escript::Data& d_dirac,const escript::Data& y_dirac) const
510 {
511 throw RipleyException("addPDEToSystem() not implemented");
512 }
513
514 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
515 const escript::Data& D, const escript::Data& d,
516 const escript::Data& d_dirac, const bool useHRZ) const
517 {
518 throw RipleyException("addPDEToLumpedSystem() not implemented");
519 }
520
521 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
522 const escript::Data& Y, const escript::Data& y,
523 const escript::Data& y_contact, const escript::Data& y_dirac) const
524 {
525 throw RipleyException("addPDEToRHS() not implemented");
526 }
527
528 void RipleyDomain::addPDEToTransportProblem(
529 escript::AbstractTransportProblem& tp,
530 escript::Data& source, const escript::Data& M,
531 const escript::Data& A, const escript::Data& B, const escript::Data& C,
532 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
533 const escript::Data& d, const escript::Data& y,
534 const escript::Data& d_contact, const escript::Data& y_contact,
535 const escript::Data& d_dirac, const escript::Data& y_dirac) const
536 {
537 throw RipleyException("addPDEToTransportProblem() not implemented");
538 }
539
540 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
541 const int blocksize, const escript::FunctionSpace& functionspace,
542 const int type) const
543 {
544 throw RipleyException("newTransportProblem() not implemented");
545 }
546
547 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
548 bool reducedColOrder) const
549 {
550 throw RipleyException("getPattern() not implemented");
551 }
552
553 dim_t RipleyDomain::getNumDataPointsGlobal() const
554 {
555 throw RipleyException("getNumDataPointsGlobal() not implemented");
556 }
557
558 IndexVector RipleyDomain::getNumNodesPerDim() const
559 {
560 throw RipleyException("getNumNodesPerDim() not implemented");
561 }
562
563 IndexVector RipleyDomain::getNumElementsPerDim() const
564 {
565 throw RipleyException("getNumElementsPerDim() not implemented");
566 }
567
568 IndexVector RipleyDomain::getNumFacesPerBoundary() const
569 {
570 throw RipleyException("getNumFacesPerBoundary() not implemented");
571 }
572
573 IndexVector RipleyDomain::getNodeDistribution() const
574 {
575 throw RipleyException("getNodeDistribution() not implemented");
576 }
577
578 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
579 {
580 throw RipleyException("getFirstCoordAndSpacing() not implemented");
581 }
582
583 dim_t RipleyDomain::getNumFaceElements() const
584 {
585 throw RipleyException("getNumFaceElements() not implemented");
586 }
587
588 dim_t RipleyDomain::getNumElements() const
589 {
590 throw RipleyException("getNumElements() not implemented");
591 }
592
593 dim_t RipleyDomain::getNumNodes() const
594 {
595 throw RipleyException("getNumNodes() not implemented");
596 }
597
598 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
599 {
600 throw RipleyException("assembleCoordinates() not implemented");
601 }
602
603 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
604 {
605 throw RipleyException("interpolateNodesOnElements() not implemented");
606 }
607
608 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
609 {
610 throw RipleyException("interpolateNodesOnFaces() not implemented");
611 }
612
613
614 } // end of namespace ripley
615

  ViewVC Help
Powered by ViewVC 1.1.26