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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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 return pair<int,int>(1, getNumNodes());
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 /*
100 case ReducedNodes:
101 return pair<int,int>(1, getNumReducedNodes());
102 case ReducedElements:
103 return pair<int,int>(1, getNumReducedElements());
104 case ReducedFaceElements:
105 return pair<int,int>(1, getNumReducedFaceElements());
106 case Points:
107 return pair<int,int>(1, getNumNodes());
108 case ReducedDegreesOfFreedom:
109 return pair<int,int>(1, getNumReducedNodes());
110 */
111 default:
112 break;
113 }
114
115 stringstream msg;
116 msg << "Invalid function space type " << fsType << " for "
117 << getDescription();
118 throw RipleyException(msg.str());
119 }
120
121 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
122 {
123 throw RipleyException("getTagFromSampleNo() not implemented");
124 }
125
126 string RipleyDomain::showTagNames() const
127 {
128 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 }
134 return ret.str();
135 }
136
137 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
138 {
139 /*
140 The idea is to use equivalence classes (i.e. types which can be
141 interpolated back and forth):
142 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 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 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 if (hasline[0]==1)
213 resultcode=Points;
214 else if (hasline[1]==1) {
215 if (hasclass[5]==1)
216 resultcode=ReducedElements;
217 else
218 resultcode=Elements;
219 } else if (hasline[2]==1) {
220 if (hasclass[7]==1)
221 resultcode=ReducedFaceElements;
222 else
223 resultcode=FaceElements;
224 } else
225 throw RipleyException("Internal Ripley Error!");
226 } else { // totlines==0
227 if (hasclass[2]==1)
228 // something from class 2
229 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
230 else // something from class 1
231 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
232 }
233 return true;
234 }
235
236 escript::Data RipleyDomain::getX() const
237 {
238 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 }
268 }
269
270 bool RipleyDomain::isCellOriented(int fsType) const
271 {
272 switch(fsType) {
273 case Nodes:
274 case DegreesOfFreedom:
275 case ReducedDegreesOfFreedom:
276 return false;
277 case Elements:
278 case FaceElements:
279 case Points:
280 case ReducedElements:
281 case ReducedFaceElements:
282 return true;
283 default:
284 break;
285 }
286 stringstream msg;
287 msg << "Illegal function space type " << fsType << " on " << getDescription();
288 throw RipleyException(msg.str());
289 }
290
291 void RipleyDomain::Print_Mesh_Info(const bool full) const
292 {
293 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
297 // 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 }
306 }
307
308 int RipleyDomain::getSystemMatrixTypeId(const int solver,
309 const int preconditioner, const int package, const bool symmetry) const
310 {
311 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
312 package, symmetry, m_mpiInfo);
313 }
314
315 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
316 const int package, const bool symmetry) const
317 {
318 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
319 package, symmetry, m_mpiInfo);
320 }
321
322 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 {
328 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 }
357
358 //
359 // the methods that follow have to be implemented by the subclasses
360 //
361
362 string RipleyDomain::getDescription() const
363 {
364 throw RipleyException("getDescription() not implemented");
365 }
366
367 bool RipleyDomain::operator==(const AbstractDomain& other) const
368 {
369 throw RipleyException("operator==() not implemented");
370 }
371
372 void RipleyDomain::write(const string& filename) const
373 {
374 throw RipleyException("write() not implemented");
375 }
376
377 void RipleyDomain::dump(const string& filename) const
378 {
379 throw RipleyException("dump() not implemented");
380 }
381
382 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
383 {
384 throw RipleyException("borrowSampleReferenceIDs() not implemented");
385 }
386
387 void RipleyDomain::setNewX(const escript::Data& arg)
388 {
389 throw RipleyException("setNewX(): Operation not supported");
390 }
391
392 void RipleyDomain::interpolateOnDomain(escript::Data& target,
393 const escript::Data& in) const
394 {
395 throw RipleyException("interpolateOnDomain() not implemented");
396 }
397
398 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
399 int fsType_target) const
400 {
401 throw RipleyException("probeInterpolationOnDomain() not implemented");
402 }
403
404 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
405 {
406 throw RipleyException("interpolateACross() not implemented");
407 }
408
409 bool RipleyDomain::probeInterpolationACross(int fsType_source,
410 const escript::AbstractDomain&, int fsType_target) const
411 {
412 throw RipleyException("probeInterpolationACross() not implemented");
413 }
414
415 void RipleyDomain::setToNormal(escript::Data& normal) const
416 {
417 throw RipleyException("setToNormal() not implemented");
418 }
419
420 void RipleyDomain::setToSize(escript::Data& size) const
421 {
422 throw RipleyException("setToSize() not implemented");
423 }
424
425 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
426 {
427 throw RipleyException("setToGradient() not implemented");
428 }
429
430 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
431 {
432 throw RipleyException("setToIntegrals() not implemented");
433 }
434
435 bool RipleyDomain::ownSample(int fsType, index_t id) const
436 {
437 throw RipleyException("ownSample() not implemented");
438 }
439
440 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const
441 {
442 throw RipleyException("setTags() not implemented");
443 }
444
445 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
446 {
447 throw RipleyException("getNumberOfTagsInUse() not implemented");
448 }
449
450 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
451 {
452 throw RipleyException("borrowListOfTagsInUse() not implemented");
453 }
454
455 bool RipleyDomain::canTag(int fsType) const
456 {
457 throw RipleyException("canTag() not implemented");
458 }
459
460
461 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
472 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 }
478
479 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 {
483 throw RipleyException("addPDEToRHS() not implemented");
484 }
485
486 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 {
495 throw RipleyException("addPDEToTransportProblem() not implemented");
496 }
497
498 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
499 const int blocksize, const escript::FunctionSpace& functionspace,
500 const int type) const
501 {
502 throw RipleyException("newTransportProblem() not implemented");
503 }
504
505 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
506 bool reducedColOrder) const
507 {
508 throw RipleyException("getPattern() not implemented");
509 }
510
511 dim_t RipleyDomain::getNumDataPointsGlobal() const
512 {
513 throw RipleyException("getNumDataPointsGlobal() not implemented");
514 }
515
516 dim_t RipleyDomain::getNumFaceElements() const
517 {
518 throw RipleyException("getNumFaceElements() not implemented");
519 }
520
521 dim_t RipleyDomain::getNumElements() const
522 {
523 throw RipleyException("getNumElements() not implemented");
524 }
525
526 dim_t RipleyDomain::getNumNodes() const
527 {
528 throw RipleyException("getNumNodes() not implemented");
529 }
530
531 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
532 {
533 throw RipleyException("assembleCoordinates() not implemented");
534 }
535
536
537 } // end of namespace ripley
538

  ViewVC Help
Powered by ViewVC 1.1.26