/[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 3702 - (show annotations)
Fri Dec 2 06:12:32 2011 UTC (7 years, 10 months ago) by caltinay
File size: 19398 byte(s)
Gradient for Rectangle elements.

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 << "getDataShape(): Unsupported function space type "
115 << functionSpaceTypeAsString(fsType) << " for " << 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 void RipleyDomain::interpolateOnDomain(escript::Data& target,
235 const escript::Data& in) const
236 {
237 const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
238 const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
239 if (inDomain != *this)
240 throw RipleyException("Illegal domain of interpolant");
241 if (targetDomain != *this)
242 throw RipleyException("Illegal domain of interpolation target");
243
244 stringstream msg;
245 msg << "interpolateOnDomain() not implemented for function space "
246 << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
247 << " -> "
248 << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
249
250 switch (in.getFunctionSpace().getTypeCode()) {
251 case Nodes:
252 case DegreesOfFreedom:
253 switch (target.getFunctionSpace().getTypeCode()) {
254 case DegreesOfFreedom:
255 target=in;
256 break;
257
258 case Elements:
259 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in));
260 break;
261
262 case FaceElements:
263 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in));
264 break;
265
266 default:
267 throw RipleyException(msg.str());
268 }
269 break;
270 default:
271 throw RipleyException(msg.str());
272 }
273 }
274
275 escript::Data RipleyDomain::getX() const
276 {
277 return escript::continuousFunction(*this).getX();
278 }
279
280 escript::Data RipleyDomain::getNormal() const
281 {
282 return escript::functionOnBoundary(*this).getNormal();
283 }
284
285 escript::Data RipleyDomain::getSize() const
286 {
287 return escript::function(*this).getSize();
288 }
289
290 void RipleyDomain::setToX(escript::Data& arg) const
291 {
292 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
293 *(arg.getFunctionSpace().getDomain()));
294 if (argDomain != *this)
295 throw RipleyException("setToX: Illegal domain of data point locations");
296 if (!arg.isExpanded())
297 throw RipleyException("setToX: Expanded Data object expected");
298
299 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
300 assembleCoordinates(arg);
301 } else {
302 // interpolate the result
303 escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
304 assembleCoordinates(contData);
305 interpolateOnDomain(arg, contData);
306 }
307 }
308
309 bool RipleyDomain::isCellOriented(int fsType) const
310 {
311 switch(fsType) {
312 case Nodes:
313 case DegreesOfFreedom:
314 case ReducedDegreesOfFreedom:
315 return false;
316 case Elements:
317 case FaceElements:
318 case Points:
319 case ReducedElements:
320 case ReducedFaceElements:
321 return true;
322 default:
323 break;
324 }
325 stringstream msg;
326 msg << "Illegal function space type " << fsType << " on " << getDescription();
327 throw RipleyException(msg.str());
328 }
329
330 void RipleyDomain::Print_Mesh_Info(const bool full) const
331 {
332 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
333 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
334 cout << "Number of dimensions: " << m_numDim << endl;
335
336 // write tags
337 if (m_tagMap.size() > 0) {
338 cout << "Tags:" << endl;
339 TagMap::const_iterator it;
340 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
341 cout << " " << setw(5) << it->second << " "
342 << it->first << endl;
343 }
344 }
345 }
346
347 int RipleyDomain::getSystemMatrixTypeId(const int solver,
348 const int preconditioner, const int package, const bool symmetry) const
349 {
350 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
351 package, symmetry, m_mpiInfo);
352 }
353
354 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
355 const int package, const bool symmetry) const
356 {
357 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
358 package, symmetry, m_mpiInfo);
359 }
360
361 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
362 const escript::FunctionSpace& row_functionspace,
363 const int column_blocksize,
364 const escript::FunctionSpace& column_functionspace,
365 const int type) const
366 {
367 bool reduceRowOrder=false;
368 bool reduceColOrder=false;
369 // is the domain right?
370 const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
371 if (row_domain!=*this)
372 throw RipleyException("Domain of row function space does not match the domain of matrix generator");
373 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
374 if (col_domain!=*this)
375 throw RipleyException("Domain of column function space does not match the domain of matrix generator");
376 // is the function space type right?
377 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
378 reduceRowOrder=true;
379 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
380 throw RipleyException("Illegal function space type for system matrix rows");
381 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
382 reduceColOrder=true;
383 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
384 throw RipleyException("Illegal function space type for system matrix columns");
385
386 // generate matrix
387 Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
388 Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
389 row_blocksize, column_blocksize, FALSE);
390 paso::checkPasoError();
391 Paso_SystemMatrixPattern_free(pattern);
392 escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
393 row_functionspace, column_blocksize, column_functionspace));
394 return sma;
395 }
396
397 void RipleyDomain::setNewX(const escript::Data& arg)
398 {
399 throw RipleyException("setNewX(): Operation not supported");
400 }
401
402 //
403 // the methods that follow have to be implemented by the subclasses
404 //
405
406 string RipleyDomain::getDescription() const
407 {
408 throw RipleyException("getDescription() not implemented");
409 }
410
411 bool RipleyDomain::operator==(const AbstractDomain& other) const
412 {
413 throw RipleyException("operator==() not implemented");
414 }
415
416 void RipleyDomain::write(const string& filename) const
417 {
418 throw RipleyException("write() not implemented");
419 }
420
421 void RipleyDomain::dump(const string& filename) const
422 {
423 throw RipleyException("dump() not implemented");
424 }
425
426 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
427 {
428 throw RipleyException("borrowSampleReferenceIDs() not implemented");
429 }
430
431 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
432 int fsType_target) const
433 {
434 throw RipleyException("probeInterpolationOnDomain() not implemented");
435 }
436
437 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
438 {
439 throw RipleyException("interpolateACross() not implemented");
440 }
441
442 bool RipleyDomain::probeInterpolationACross(int fsType_source,
443 const escript::AbstractDomain&, int fsType_target) const
444 {
445 throw RipleyException("probeInterpolationACross() not implemented");
446 }
447
448 void RipleyDomain::setToNormal(escript::Data& normal) const
449 {
450 throw RipleyException("setToNormal() not implemented");
451 }
452
453 void RipleyDomain::setToSize(escript::Data& size) const
454 {
455 throw RipleyException("setToSize() not implemented");
456 }
457
458 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
459 {
460 throw RipleyException("setToGradient() not implemented");
461 }
462
463 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
464 {
465 throw RipleyException("setToIntegrals() not implemented");
466 }
467
468 bool RipleyDomain::ownSample(int fsType, index_t id) const
469 {
470 throw RipleyException("ownSample() not implemented");
471 }
472
473 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const
474 {
475 throw RipleyException("setTags() not implemented");
476 }
477
478 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
479 {
480 throw RipleyException("getNumberOfTagsInUse() not implemented");
481 }
482
483 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
484 {
485 throw RipleyException("borrowListOfTagsInUse() not implemented");
486 }
487
488 bool RipleyDomain::canTag(int fsType) const
489 {
490 throw RipleyException("canTag() not implemented");
491 }
492
493
494 void RipleyDomain::addPDEToSystem(
495 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
496 const escript::Data& A, const escript::Data& B, const escript::Data& C,
497 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
498 const escript::Data& d, const escript::Data& y,
499 const escript::Data& d_contact, const escript::Data& y_contact,
500 const escript::Data& d_dirac,const escript::Data& y_dirac) const
501 {
502 throw RipleyException("addPDEToSystem() not implemented");
503 }
504
505 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
506 const escript::Data& D, const escript::Data& d,
507 const escript::Data& d_dirac, const bool useHRZ) const
508 {
509 throw RipleyException("addPDEToLumpedSystem() not implemented");
510 }
511
512 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
513 const escript::Data& Y, const escript::Data& y,
514 const escript::Data& y_contact, const escript::Data& y_dirac) const
515 {
516 throw RipleyException("addPDEToRHS() not implemented");
517 }
518
519 void RipleyDomain::addPDEToTransportProblem(
520 escript::AbstractTransportProblem& tp,
521 escript::Data& source, const escript::Data& M,
522 const escript::Data& A, const escript::Data& B, const escript::Data& C,
523 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
524 const escript::Data& d, const escript::Data& y,
525 const escript::Data& d_contact, const escript::Data& y_contact,
526 const escript::Data& d_dirac, const escript::Data& y_dirac) const
527 {
528 throw RipleyException("addPDEToTransportProblem() not implemented");
529 }
530
531 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
532 const int blocksize, const escript::FunctionSpace& functionspace,
533 const int type) const
534 {
535 throw RipleyException("newTransportProblem() not implemented");
536 }
537
538 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
539 bool reducedColOrder) const
540 {
541 throw RipleyException("getPattern() not implemented");
542 }
543
544 dim_t RipleyDomain::getNumDataPointsGlobal() const
545 {
546 throw RipleyException("getNumDataPointsGlobal() not implemented");
547 }
548
549 IndexVector RipleyDomain::getNumNodesPerDim() const
550 {
551 throw RipleyException("getNumNodesPerDim() not implemented");
552 }
553
554 IndexVector RipleyDomain::getNumElementsPerDim() const
555 {
556 throw RipleyException("getNumElementsPerDim() not implemented");
557 }
558
559 IndexVector RipleyDomain::getNumFacesPerBoundary() const
560 {
561 throw RipleyException("getNumFacesPerBoundary() not implemented");
562 }
563
564 IndexVector RipleyDomain::getNodeDistribution() const
565 {
566 throw RipleyException("getNodeDistribution() not implemented");
567 }
568
569 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
570 {
571 throw RipleyException("getFirstCoordAndSpacing() not implemented");
572 }
573
574 dim_t RipleyDomain::getNumFaceElements() const
575 {
576 throw RipleyException("getNumFaceElements() not implemented");
577 }
578
579 dim_t RipleyDomain::getNumElements() const
580 {
581 throw RipleyException("getNumElements() not implemented");
582 }
583
584 dim_t RipleyDomain::getNumNodes() const
585 {
586 throw RipleyException("getNumNodes() not implemented");
587 }
588
589 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
590 {
591 throw RipleyException("assembleCoordinates() not implemented");
592 }
593
594 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in) const
595 {
596 throw RipleyException("interpolateNodesOnElements() not implemented");
597 }
598
599 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in) const
600 {
601 throw RipleyException("interpolateNodesOnFaces() not implemented");
602 }
603
604
605 } // end of namespace ripley
606

  ViewVC Help
Powered by ViewVC 1.1.26