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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3766 - (show annotations)
Thu Jan 12 08:17:49 2012 UTC (7 years, 8 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 45120 byte(s)
Checkpoint. Getting weipa to like ripley again.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2012 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::operator==(const AbstractDomain& other) const
51 {
52 const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
53 if (o) {
54 return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
55 && m_elementTags==o->m_elementTags
56 && m_faceTags==o->m_faceTags);
57 }
58 return false;
59 }
60
61 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
62 {
63 switch (fsType) {
64 case DegreesOfFreedom:
65 case ReducedDegreesOfFreedom:
66 case Nodes:
67 case ReducedNodes:
68 case Elements:
69 case ReducedElements:
70 case FaceElements:
71 case ReducedFaceElements:
72 case Points:
73 return true;
74 default:
75 break;
76 }
77 return false;
78 }
79
80 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
81 {
82 switch (fsType) {
83 case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84 case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85 case Nodes: return "Ripley_Nodes";
86 case ReducedNodes: return "Ripley_ReducedNodes";
87 case Elements: return "Ripley_Elements";
88 case ReducedElements: return "Ripley_ReducedElements";
89 case FaceElements: return "Ripley_FaceElements";
90 case ReducedFaceElements: return "Ripley_ReducedFaceElements";
91 case Points: return "Ripley_Points";
92 default:
93 break;
94 }
95 return "Invalid function space type code";
96 }
97
98 pair<int,int> RipleyDomain::getDataShape(int fsType) const
99 {
100 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
101 switch (fsType) {
102 case Nodes:
103 return pair<int,int>(1, getNumNodes());
104 case ReducedNodes: //FIXME: reduced
105 if (getCoarseMesh())
106 return getCoarseMesh()->getDataShape(Nodes);
107 break;
108 case DegreesOfFreedom:
109 case ReducedDegreesOfFreedom: //FIXME: reduced
110 return pair<int,int>(1, getNumDOF());
111 case Elements:
112 return pair<int,int>(ptsPerSample, getNumElements());
113 case FaceElements:
114 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
115 case ReducedElements:
116 return pair<int,int>(1, getNumElements());
117 case ReducedFaceElements:
118 return pair<int,int>(1, getNumFaceElements());
119 case Points:
120 return pair<int,int>(1, 0); //FIXME: dirac
121 default:
122 break;
123 }
124
125 stringstream msg;
126 msg << "getDataShape(): Unsupported function space type " << fsType
127 << " for " << getDescription();
128 throw RipleyException(msg.str());
129 }
130
131 string RipleyDomain::showTagNames() const
132 {
133 stringstream ret;
134 TagMap::const_iterator it;
135 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
136 if (it!=m_tagMap.begin()) ret << ", ";
137 ret << it->first;
138 }
139 return ret.str();
140 }
141
142 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
143 {
144 /*
145 The idea is to use equivalence classes (i.e. types which can be
146 interpolated back and forth):
147 class 0: DOF <-> Nodes
148 class 1: ReducedDOF <-> ReducedNodes
149 class 2: Points
150 class 3: Elements
151 class 4: ReducedElements
152 class 5: FaceElements
153 class 6: ReducedFaceElements
154
155 There is also a set of lines. Interpolation is possible down a line but not
156 between lines.
157 class 0 and 1 belong to all lines so aren't considered.
158 line 0: class 2
159 line 1: class 3,4
160 line 2: class 5,6
161
162 For classes with multiple members (eg class 1) we have vars to record if
163 there is at least one instance. e.g. hasnodes is true if we have at least
164 one instance of Nodes.
165 */
166 if (fs.empty())
167 return false;
168 vector<bool> hasclass(7, false);
169 vector<int> hasline(3, 0);
170 bool hasnodes=false;
171 bool hasrednodes=false;
172 for (size_t i=0; i<fs.size(); ++i) {
173 switch (fs[i]) {
174 case Nodes: hasnodes=true; // fall through
175 case DegreesOfFreedom:
176 hasclass[0]=true;
177 break;
178 case ReducedNodes: hasrednodes=true; // fall through
179 case ReducedDegreesOfFreedom:
180 hasclass[1]=true;
181 break;
182 case Points:
183 hasline[0]=1;
184 hasclass[2]=true;
185 break;
186 case Elements:
187 hasclass[3]=true;
188 hasline[1]=1;
189 break;
190 case ReducedElements:
191 hasclass[4]=true;
192 hasline[1]=1;
193 break;
194 case FaceElements:
195 hasclass[5]=true;
196 hasline[2]=1;
197 break;
198 case ReducedFaceElements:
199 hasclass[6]=true;
200 hasline[2]=1;
201 break;
202 default:
203 return false;
204 }
205 }
206 int numLines=hasline[0]+hasline[1]+hasline[2];
207
208 // fail if we have more than one leaf group
209 // = there are at least two branches we can't interpolate between
210 if (numLines > 1)
211 return false;
212 else if (numLines==1) {
213 // we have points
214 if (hasline[0]==1)
215 resultcode=Points;
216 else if (hasline[1]==1) {
217 if (hasclass[4])
218 resultcode=ReducedElements;
219 else
220 resultcode=Elements;
221 } else { // hasline[2]==1
222 if (hasclass[6])
223 resultcode=ReducedFaceElements;
224 else
225 resultcode=FaceElements;
226 }
227 } else { // numLines==0
228 if (hasclass[1])
229 // something from class 1
230 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
231 else // something from class 0
232 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
233 }
234 return true;
235 }
236
237 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
238 int fsType_target) const
239 {
240 if (!isValidFunctionSpaceType(fsType_target)) {
241 stringstream msg;
242 msg << "probeInterpolationOnDomain(): Invalid functionspace type "
243 << fsType_target << " for " << getDescription();
244 throw RipleyException(msg.str());
245 }
246
247 switch (fsType_source) {
248 case Nodes:
249 case DegreesOfFreedom:
250 return true;
251 case ReducedNodes:
252 case ReducedDegreesOfFreedom:
253 return (fsType_target != Nodes &&
254 fsType_target != DegreesOfFreedom);
255 case Elements:
256 return (fsType_target==Elements ||
257 fsType_target==ReducedElements);
258 case FaceElements:
259 return (fsType_target==FaceElements ||
260 fsType_target==ReducedFaceElements);
261 case ReducedElements:
262 case ReducedFaceElements:
263 case Points:
264 return (fsType_target==fsType_source);
265
266 default: {
267 stringstream msg;
268 msg << "probeInterpolationOnDomain(): Invalid functionspace type "
269 << fsType_source << " for " << getDescription();
270 throw RipleyException(msg.str());
271 }
272 }
273 }
274
275 void RipleyDomain::interpolateOnDomain(escript::Data& target,
276 const escript::Data& in) const
277 {
278 const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
279 const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
280 if (inDomain != *this)
281 throw RipleyException("Illegal domain of interpolant");
282 if (targetDomain != *this)
283 throw RipleyException("Illegal domain of interpolation target");
284
285 stringstream msg;
286 msg << "interpolateOnDomain() not implemented for function space "
287 << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
288 << " -> "
289 << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
290
291 const int inFS = in.getFunctionSpace().getTypeCode();
292 const int outFS = target.getFunctionSpace().getTypeCode();
293
294 // simplest case: 1:1 copy
295 if (inFS==outFS) {
296 copyData(target, *const_cast<escript::Data*>(&in));
297 // not allowed: reduced->non-reduced
298 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
299 && (outFS==Nodes || outFS==DegreesOfFreedom)) {
300 throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
301 } else if ((inFS==Elements && outFS==ReducedElements)
302 || (inFS==FaceElements && outFS==ReducedFaceElements)) {
303 averageData(target, *const_cast<escript::Data*>(&in));
304 } else {
305 switch (inFS) {
306 case Nodes:
307 case ReducedNodes: //FIXME: reduced
308 switch (outFS) {
309 case DegreesOfFreedom:
310 case ReducedDegreesOfFreedom: //FIXME: reduced
311 if (getMPISize()==1)
312 copyData(target, *const_cast<escript::Data*>(&in));
313 else
314 nodesToDOF(target,*const_cast<escript::Data*>(&in));
315 break;
316
317 case Nodes:
318 copyData(target, *const_cast<escript::Data*>(&in));
319 case ReducedNodes: //FIXME: reduced
320 break;
321
322 case Elements:
323 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
324 break;
325
326 case ReducedElements:
327 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
328 break;
329
330 case FaceElements:
331 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
332 break;
333
334 case ReducedFaceElements:
335 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
336 break;
337
338 default:
339 throw RipleyException(msg.str());
340 }
341 break;
342
343 case DegreesOfFreedom:
344 case ReducedDegreesOfFreedom: //FIXME: reduced
345 switch (outFS) {
346 case Nodes:
347 case ReducedNodes: //FIXME: reduced
348 if (getMPISize()==1)
349 copyData(target, *const_cast<escript::Data*>(&in));
350 else
351 dofToNodes(target, *const_cast<escript::Data*>(&in));
352 break;
353
354 case DegreesOfFreedom:
355 case ReducedDegreesOfFreedom: //FIXME: reduced
356 copyData(target, *const_cast<escript::Data*>(&in));
357 break;
358
359 case Elements:
360 case ReducedElements:
361 if (getMPISize()==1) {
362 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
363 } else {
364 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
365 escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
366 interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
367 }
368 break;
369
370 case FaceElements:
371 case ReducedFaceElements:
372 if (getMPISize()==1) {
373 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
374 } else {
375 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
376 escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
377 interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
378 }
379 break;
380
381 default:
382 throw RipleyException(msg.str());
383 }
384 break;
385 default:
386 throw RipleyException(msg.str());
387 }
388 }
389 }
390
391 escript::Data RipleyDomain::getX() const
392 {
393 return escript::continuousFunction(*this).getX();
394 }
395
396 escript::Data RipleyDomain::getNormal() const
397 {
398 return escript::functionOnBoundary(*this).getNormal();
399 }
400
401 escript::Data RipleyDomain::getSize() const
402 {
403 return escript::function(*this).getSize();
404 }
405
406 void RipleyDomain::setToX(escript::Data& arg) const
407 {
408 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
409 *(arg.getFunctionSpace().getDomain()));
410 if (argDomain != *this)
411 throw RipleyException("setToX: Illegal domain of data point locations");
412 if (!arg.isExpanded())
413 throw RipleyException("setToX: Expanded Data object expected");
414
415 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
416 assembleCoordinates(arg);
417 } else {
418 // interpolate the result
419 escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
420 assembleCoordinates(contData);
421 interpolateOnDomain(arg, contData);
422 }
423 }
424
425 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
426 {
427 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
428 *(arg.getFunctionSpace().getDomain()));
429 if (argDomain != *this)
430 throw RipleyException("setToGradient: Illegal domain of gradient argument");
431 const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
432 *(grad.getFunctionSpace().getDomain()));
433 if (gradDomain != *this)
434 throw RipleyException("setToGradient: Illegal domain of gradient");
435
436 switch (grad.getFunctionSpace().getTypeCode()) {
437 case Elements:
438 case ReducedElements:
439 case FaceElements:
440 case ReducedFaceElements:
441 break;
442 default: {
443 stringstream msg;
444 msg << "setToGradient(): not supported for "
445 << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
446 throw RipleyException(msg.str());
447 }
448 }
449
450 if (getMPISize()>1) {
451 if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
452 escript::Data contArg(arg, escript::continuousFunction(*this));
453 assembleGradient(grad, contArg);
454 } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
455 escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
456 assembleGradient(grad, contArg);
457 } else {
458 assembleGradient(grad, *const_cast<escript::Data*>(&arg));
459 }
460 } else {
461 assembleGradient(grad, *const_cast<escript::Data*>(&arg));
462 }
463 }
464
465 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
466 {
467 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
468 *(arg.getFunctionSpace().getDomain()));
469 if (argDomain != *this)
470 throw RipleyException("setToIntegrals: Illegal domain of integration kernel");
471
472 switch (arg.getFunctionSpace().getTypeCode()) {
473 case Nodes:
474 case ReducedNodes:
475 case DegreesOfFreedom:
476 case ReducedDegreesOfFreedom:
477 {
478 escript::Data funcArg(arg, escript::function(*this));
479 assembleIntegrate(integrals, funcArg);
480 }
481 break;
482 case Elements:
483 case ReducedElements:
484 case FaceElements:
485 case ReducedFaceElements:
486 assembleIntegrate(integrals, *const_cast<escript::Data*>(&arg));
487 break;
488 default: {
489 stringstream msg;
490 msg << "setToIntegrals(): not supported for "
491 << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
492 throw RipleyException(msg.str());
493 }
494 }
495
496 }
497
498 bool RipleyDomain::isCellOriented(int fsType) const
499 {
500 switch(fsType) {
501 case Nodes:
502 case ReducedNodes:
503 case DegreesOfFreedom:
504 case ReducedDegreesOfFreedom:
505 return false;
506 case Elements:
507 case FaceElements:
508 case Points:
509 case ReducedElements:
510 case ReducedFaceElements:
511 return true;
512 default:
513 break;
514 }
515 stringstream msg;
516 msg << "isCellOriented(): Illegal function space type " << fsType
517 << " on " << getDescription();
518 throw RipleyException(msg.str());
519 }
520
521 bool RipleyDomain::canTag(int fsType) const
522 {
523 switch(fsType) {
524 case Nodes:
525 case Elements:
526 case ReducedElements:
527 case FaceElements:
528 case ReducedFaceElements:
529 return true;
530 case DegreesOfFreedom:
531 case ReducedDegreesOfFreedom:
532 case Points:
533 case ReducedNodes:
534 return false;
535 default:
536 break;
537 }
538 stringstream msg;
539 msg << "canTag(): Illegal function space type " << fsType << " on "
540 << getDescription();
541 throw RipleyException(msg.str());
542 }
543
544 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
545 {
546 IndexVector* target=NULL;
547 dim_t num=0;
548
549 switch(fsType) {
550 case Nodes:
551 num=getNumNodes();
552 target=&m_nodeTags;
553 break;
554 case Elements:
555 case ReducedElements:
556 num=getNumElements();
557 target=&m_elementTags;
558 break;
559 case FaceElements:
560 case ReducedFaceElements:
561 num=getNumFaceElements();
562 target=&m_faceTags;
563 break;
564 default: {
565 stringstream msg;
566 msg << "setTags(): not implemented for "
567 << functionSpaceTypeAsString(fsType);
568 throw RipleyException(msg.str());
569 }
570 }
571 if (target->size() != num) {
572 target->assign(num, -1);
573 }
574
575 escript::Data& mask=*const_cast<escript::Data*>(&cMask);
576 #pragma omp parallel for
577 for (index_t i=0; i<num; i++) {
578 if (mask.getSampleDataRO(i)[0] > 0) {
579 (*target)[i]=newTag;
580 }
581 }
582 updateTagsInUse(fsType);
583 }
584
585 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
586 {
587 switch(fsType) {
588 case Nodes:
589 if (m_nodeTags.size() > sampleNo)
590 return m_nodeTags[sampleNo];
591 break;
592 case Elements:
593 case ReducedElements:
594 if (m_elementTags.size() > sampleNo)
595 return m_elementTags[sampleNo];
596 break;
597 case FaceElements:
598 case ReducedFaceElements:
599 if (m_faceTags.size() > sampleNo)
600 return m_faceTags[sampleNo];
601 break;
602 default: {
603 stringstream msg;
604 msg << "getTagFromSampleNo(): not implemented for "
605 << functionSpaceTypeAsString(fsType);
606 throw RipleyException(msg.str());
607 }
608 }
609 return -1;
610 }
611
612 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
613 {
614 switch(fsType) {
615 case Nodes:
616 return m_nodeTagsInUse.size();
617 case Elements:
618 case ReducedElements:
619 return m_elementTagsInUse.size();
620 case FaceElements:
621 case ReducedFaceElements:
622 return m_faceTagsInUse.size();
623 default: {
624 stringstream msg;
625 msg << "getNumberOfTagsInUse(): not implemented for "
626 << functionSpaceTypeAsString(fsType);
627 throw RipleyException(msg.str());
628 }
629 }
630 }
631
632 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
633 {
634 switch(fsType) {
635 case Nodes:
636 return &m_nodeTagsInUse[0];
637 case Elements:
638 case ReducedElements:
639 return &m_elementTagsInUse[0];
640 case FaceElements:
641 case ReducedFaceElements:
642 return &m_faceTagsInUse[0];
643 default: {
644 stringstream msg;
645 msg << "borrowListOfTagsInUse(): not implemented for "
646 << functionSpaceTypeAsString(fsType);
647 throw RipleyException(msg.str());
648 }
649 }
650 }
651
652 void RipleyDomain::Print_Mesh_Info(const bool full) const
653 {
654 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
655 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
656 cout << "Number of dimensions: " << m_numDim << endl;
657
658 // write tags
659 if (m_tagMap.size() > 0) {
660 cout << "Tags:" << endl;
661 TagMap::const_iterator it;
662 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
663 cout << " " << setw(5) << it->second << " "
664 << it->first << endl;
665 }
666 }
667 }
668
669 int RipleyDomain::getSystemMatrixTypeId(const int solver,
670 const int preconditioner, const int package, const bool symmetry) const
671 {
672 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
673 package, symmetry, m_mpiInfo);
674 }
675
676 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
677 const int package, const bool symmetry) const
678 {
679 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
680 package, symmetry, m_mpiInfo);
681 }
682
683 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
684 const escript::FunctionSpace& row_functionspace,
685 const int column_blocksize,
686 const escript::FunctionSpace& column_functionspace,
687 const int type) const
688 {
689 bool reduceRowOrder=false;
690 bool reduceColOrder=false;
691 // is the domain right?
692 const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
693 if (row_domain != *this)
694 throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
695 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
696 if (col_domain != *this)
697 throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
698 // is the function space type right?
699 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
700 reduceRowOrder=true;
701 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
702 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
703 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
704 reduceColOrder=true;
705 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
706 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
707
708 // generate matrix
709 Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
710 Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
711 row_blocksize, column_blocksize, FALSE);
712 paso::checkPasoError();
713 Paso_SystemMatrixPattern_free(pattern);
714 escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
715 row_functionspace, column_blocksize, column_functionspace));
716 return sma;
717 }
718
719 void RipleyDomain::addPDEToSystem(
720 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
721 const escript::Data& A, const escript::Data& B, const escript::Data& C,
722 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
723 const escript::Data& d, const escript::Data& y,
724 const escript::Data& d_contact, const escript::Data& y_contact,
725 const escript::Data& d_dirac,const escript::Data& y_dirac) const
726 {
727 if (!d_contact.isEmpty() || !y_contact.isEmpty())
728 throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
729
730 paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
731 if (!sma)
732 throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
733
734 if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
735 throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
736
737 int fsType=UNKNOWN;
738 fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
739 fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
740 fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
741 fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
742 fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
743 fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
744 fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
745 fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
746 fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
747 fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
748 if (fsType==UNKNOWN)
749 return;
750
751 const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
752
753 //TODO: more input param checks (shape, etc)
754
755 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
756
757 if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
758 throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
759
760 const int numEq=S->logical_row_block_size;
761 const int numComp=S->logical_col_block_size;
762
763 if (numEq != numComp)
764 throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
765 //TODO: more system matrix checks
766
767 if (numEq==1)
768 if (reducedOrder)
769 assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y, d, y);
770 else
771 assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
772 else
773 if (reducedOrder)
774 assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y, d, y);
775 else
776 assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
777 }
778
779 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
780 const escript::Data& Y, const escript::Data& y,
781 const escript::Data& y_contact, const escript::Data& y_dirac) const
782 {
783 if (!y_contact.isEmpty())
784 throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
785
786 if (rhs.isEmpty()) {
787 if (!X.isEmpty() || !Y.isEmpty())
788 throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
789 else
790 return;
791 }
792
793 if (rhs.getDataPointSize() == 1)
794 assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
795 else
796 assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
797 }
798
799 void RipleyDomain::setNewX(const escript::Data& arg)
800 {
801 throw RipleyException("setNewX(): Operation not supported");
802 }
803
804 //protected
805 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
806 {
807 const dim_t numComp = in.getDataPointSize();
808 const dim_t dpp = in.getNumDataPointsPerSample();
809 out.requireWrite();
810 #pragma omp parallel for
811 for (index_t i=0; i<in.getNumSamples(); i++) {
812 const double* src = in.getSampleDataRO(i);
813 double* dest = out.getSampleDataRW(i);
814 for (index_t c=0; c<numComp; c++) {
815 double res=0.;
816 for (index_t q=0; q<dpp; q++)
817 res+=src[c+q*numComp];
818 *dest++ = res/dpp;
819 }
820 }
821 }
822
823 //protected
824 void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
825 {
826 const dim_t numComp = in.getDataPointSize();
827 out.requireWrite();
828 #pragma omp parallel for
829 for (index_t i=0; i<in.getNumSamples(); i++) {
830 const double* src = in.getSampleDataRO(i);
831 copy(src, src+numComp, out.getSampleDataRW(i));
832 }
833 }
834
835 //protected
836 void RipleyDomain::updateTagsInUse(int fsType) const
837 {
838 IndexVector* tagsInUse=NULL;
839 const IndexVector* tags=NULL;
840 switch(fsType) {
841 case Nodes:
842 tags=&m_nodeTags;
843 tagsInUse=&m_nodeTagsInUse;
844 break;
845 case Elements:
846 case ReducedElements:
847 tags=&m_elementTags;
848 tagsInUse=&m_elementTagsInUse;
849 break;
850 case FaceElements:
851 case ReducedFaceElements:
852 tags=&m_faceTags;
853 tagsInUse=&m_faceTagsInUse;
854 break;
855 default:
856 return;
857 }
858
859 // gather global unique tag values from tags into tagsInUse
860 tagsInUse->clear();
861 index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
862
863 while (true) {
864 // find smallest value bigger than lastFoundValue
865 minFoundValue = INDEX_T_MAX;
866 #pragma omp parallel private(local_minFoundValue)
867 {
868 local_minFoundValue = minFoundValue;
869 #pragma omp for schedule(static) nowait
870 for (size_t i = 0; i < tags->size(); i++) {
871 const index_t v = (*tags)[i];
872 if ((v > lastFoundValue) && (v < local_minFoundValue))
873 local_minFoundValue = v;
874 }
875 #pragma omp critical
876 {
877 if (local_minFoundValue < minFoundValue)
878 minFoundValue = local_minFoundValue;
879 }
880 }
881 #ifdef ESYS_MPI
882 local_minFoundValue = minFoundValue;
883 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
884 #endif
885
886 // if we found a new value add it to the tagsInUse vector
887 if (minFoundValue < INDEX_T_MAX) {
888 tagsInUse->push_back(minFoundValue);
889 lastFoundValue = minFoundValue;
890 } else
891 break;
892 }
893 }
894
895 //protected
896 Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
897 const IndexVector& index, const dim_t M, const dim_t N) const
898 {
899 // paso will manage the memory
900 index_t* indexC = MEMALLOC(index.size(), index_t);
901 index_t* ptrC = MEMALLOC(ptr.size(), index_t);
902 copy(index.begin(), index.end(), indexC);
903 copy(ptr.begin(), ptr.end(), ptrC);
904 return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
905 }
906
907 //protected
908 Paso_Pattern* RipleyDomain::createMainPattern() const
909 {
910 IndexVector ptr(1,0);
911 IndexVector index;
912
913 for (index_t i=0; i<getNumDOF(); i++) {
914 // add the DOF itself
915 index.push_back(i);
916 const dim_t num=insertNeighbourNodes(index, i);
917 ptr.push_back(ptr.back()+num+1);
918 }
919
920 return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
921 }
922
923 //protected
924 void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
925 const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
926 {
927 IndexVector ptr(1,0);
928 IndexVector index;
929 for (index_t i=0; i<getNumDOF(); i++) {
930 index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
931 ptr.push_back(ptr.back()+colIndices[i].size());
932 }
933
934 const dim_t M=ptr.size()-1;
935 *colPattern=createPasoPattern(ptr, index, M, N);
936
937 IndexVector rowPtr(1,0);
938 IndexVector rowIndex;
939 for (dim_t id=0; id<N; id++) {
940 dim_t n=0;
941 for (dim_t i=0; i<M; i++) {
942 for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
943 if (index[j]==id) {
944 rowIndex.push_back(i);
945 n++;
946 break;
947 }
948 }
949 }
950 rowPtr.push_back(rowPtr.back()+n);
951 }
952
953 // M and N are now reversed
954 *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
955 }
956
957 //protected
958 void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
959 const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
960 dim_t num_Sol, const vector<double>& array) const
961 {
962 const dim_t numMyCols = mat->pattern->mainPattern->numInput;
963 const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
964 const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
965 const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
966
967 const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
968 const index_t* mainBlock_index = mat->mainBlock->pattern->index;
969 double* mainBlock_val = mat->mainBlock->val;
970 const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
971 const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
972 double* col_coupleBlock_val = mat->col_coupleBlock->val;
973 const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
974 const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
975 double* row_coupleBlock_val = mat->row_coupleBlock->val;
976
977 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
978 // down columns of array
979 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
980 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
981 // only look at the matrix rows stored on this processor
982 if (i_row < numMyRows) {
983 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
984 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
985 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
986 if (i_col < numMyCols) {
987 for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
988 if (mainBlock_index[k] == i_col) {
989 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
990 const dim_t i_Sol=ic+mat->col_block_size*l_col;
991 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
992 const dim_t i_Eq=ir+mat->row_block_size*l_row;
993 mainBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
994 }
995 }
996 break;
997 }
998 }
999 } else {
1000 for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1001 if (col_coupleBlock_index[k] == i_col - numMyCols) {
1002 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1003 const dim_t i_Sol=ic+mat->col_block_size*l_col;
1004 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1005 const dim_t i_Eq=ir+mat->row_block_size*l_row;
1006 col_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1007 }
1008 }
1009 break;
1010 }
1011 }
1012 }
1013 }
1014 }
1015 } else {
1016 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1017 // across rows of array
1018 for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1019 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1020 if (i_col < numMyCols) {
1021 for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1022 k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1023 {
1024 if (row_coupleBlock_index[k] == i_col) {
1025 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1026 const dim_t i_Sol=ic+mat->col_block_size*l_col;
1027 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1028 const dim_t i_Eq=ir+mat->row_block_size*l_row;
1029 row_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1030 }
1031 }
1032 break;
1033 }
1034 }
1035 }
1036 }
1037 }
1038 }
1039 }
1040 }
1041 }
1042
1043 //
1044 // the methods that follow have to be implemented by the subclasses
1045 //
1046
1047 string RipleyDomain::getDescription() const
1048 {
1049 throw RipleyException("getDescription() not implemented");
1050 }
1051
1052 void RipleyDomain::write(const string& filename) const
1053 {
1054 throw RipleyException("write() not implemented");
1055 }
1056
1057 void RipleyDomain::dump(const string& filename) const
1058 {
1059 throw RipleyException("dump() not implemented");
1060 }
1061
1062 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
1063 {
1064 throw RipleyException("borrowSampleReferenceIDs() not implemented");
1065 }
1066
1067 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1068 {
1069 throw RipleyException("interpolateACross() not implemented");
1070 }
1071
1072 bool RipleyDomain::probeInterpolationACross(int fsType_source,
1073 const escript::AbstractDomain&, int fsType_target) const
1074 {
1075 throw RipleyException("probeInterpolationACross() not implemented");
1076 }
1077
1078 void RipleyDomain::setToNormal(escript::Data& normal) const
1079 {
1080 throw RipleyException("setToNormal() not implemented");
1081 }
1082
1083 void RipleyDomain::setToSize(escript::Data& size) const
1084 {
1085 throw RipleyException("setToSize() not implemented");
1086 }
1087
1088 bool RipleyDomain::ownSample(int fsType, index_t id) const
1089 {
1090 throw RipleyException("ownSample() not implemented");
1091 }
1092
1093 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1094 const escript::Data& D, const escript::Data& d,
1095 const escript::Data& d_dirac, const bool useHRZ) const
1096 {
1097 throw RipleyException("addPDEToLumpedSystem() not implemented");
1098 }
1099
1100 void RipleyDomain::addPDEToTransportProblem(
1101 escript::AbstractTransportProblem& tp,
1102 escript::Data& source, const escript::Data& M,
1103 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1104 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1105 const escript::Data& d, const escript::Data& y,
1106 const escript::Data& d_contact, const escript::Data& y_contact,
1107 const escript::Data& d_dirac, const escript::Data& y_dirac) const
1108 {
1109 throw RipleyException("addPDEToTransportProblem() not implemented");
1110 }
1111
1112 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
1113 const int blocksize, const escript::FunctionSpace& functionspace,
1114 const int type) const
1115 {
1116 throw RipleyException("newTransportProblem() not implemented");
1117 }
1118
1119 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1120 bool reducedColOrder) const
1121 {
1122 throw RipleyException("getPattern() not implemented");
1123 }
1124
1125 dim_t RipleyDomain::getNumDataPointsGlobal() const
1126 {
1127 throw RipleyException("getNumDataPointsGlobal() not implemented");
1128 }
1129
1130 IndexVector RipleyDomain::getNumNodesPerDim() const
1131 {
1132 throw RipleyException("getNumNodesPerDim() not implemented");
1133 }
1134
1135 IndexVector RipleyDomain::getNumElementsPerDim() const
1136 {
1137 throw RipleyException("getNumElementsPerDim() not implemented");
1138 }
1139
1140 IndexVector RipleyDomain::getNumFacesPerBoundary() const
1141 {
1142 throw RipleyException("getNumFacesPerBoundary() not implemented");
1143 }
1144
1145 IndexVector RipleyDomain::getNodeDistribution() const
1146 {
1147 throw RipleyException("getNodeDistribution() not implemented");
1148 }
1149
1150 IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1151 {
1152 throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1153 }
1154
1155 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1156 {
1157 throw RipleyException("getFirstCoordAndSpacing() not implemented");
1158 }
1159
1160 dim_t RipleyDomain::getNumFaceElements() const
1161 {
1162 throw RipleyException("getNumFaceElements() not implemented");
1163 }
1164
1165 dim_t RipleyDomain::getNumElements() const
1166 {
1167 throw RipleyException("getNumElements() not implemented");
1168 }
1169
1170 dim_t RipleyDomain::getNumNodes() const
1171 {
1172 throw RipleyException("getNumNodes() not implemented");
1173 }
1174
1175 dim_t RipleyDomain::getNumDOF() const
1176 {
1177 throw RipleyException("getNumDOF() not implemented");
1178 }
1179
1180 escript::Domain_ptr RipleyDomain::getCoarseMesh() const
1181 {
1182 throw RipleyException("getCoarseMesh() not implemented");
1183 }
1184
1185 dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1186 {
1187 throw RipleyException("insertNeighbourNodes() not implemented");
1188 }
1189
1190 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1191 {
1192 throw RipleyException("assembleCoordinates() not implemented");
1193 }
1194
1195 void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1196 {
1197 throw RipleyException("assembleGradient() not implemented");
1198 }
1199
1200 void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1201 {
1202 throw RipleyException("assembleIntegrate() not implemented");
1203 }
1204
1205 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1206 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1207 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1208 const escript::Data& d, const escript::Data& y) const
1209 {
1210 throw RipleyException("assemblePDESingle() not implemented");
1211 }
1212
1213 void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1214 escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1215 const escript::Data& C, const escript::Data& D, const escript::Data& X,
1216 const escript::Data& Y, const escript::Data& d, const escript::Data& y) const
1217 {
1218 throw RipleyException("assemblePDESingleReduced() not implemented");
1219 }
1220
1221 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1222 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1223 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1224 const escript::Data& d, const escript::Data& y) const
1225 {
1226 throw RipleyException("assemblePDESystem() not implemented");
1227 }
1228
1229 void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1230 escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1231 const escript::Data& C, const escript::Data& D, const escript::Data& X,
1232 const escript::Data& Y, const escript::Data& d, const escript::Data& y) const
1233 {
1234 throw RipleyException("assemblePDESystemReduced() not implemented");
1235 }
1236
1237 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1238 {
1239 throw RipleyException("interpolateNodesOnElements() not implemented");
1240 }
1241
1242 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1243 {
1244 throw RipleyException("interpolateNodesOnFaces() not implemented");
1245 }
1246
1247 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1248 {
1249 throw RipleyException("nodesToDOF() not implemented");
1250 }
1251
1252 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1253 {
1254 throw RipleyException("dofToNodes() not implemented");
1255 }
1256
1257 } // end of namespace ripley
1258

  ViewVC Help
Powered by ViewVC 1.1.26