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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3793 - (show annotations)
Wed Feb 1 07:39:43 2012 UTC (7 years, 6 months ago) by gross
File size: 45597 byte(s)
new implementation of FCT solver with some modifications to the python interface
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 case ReducedNodes: //FIXME: reduced
104 return pair<int,int>(1, getNumNodes());
105 case DegreesOfFreedom:
106 case ReducedDegreesOfFreedom: //FIXME: reduced
107 return pair<int,int>(1, getNumDOF());
108 case Elements:
109 return pair<int,int>(ptsPerSample, getNumElements());
110 case FaceElements:
111 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
112 case ReducedElements:
113 return pair<int,int>(1, getNumElements());
114 case ReducedFaceElements:
115 return pair<int,int>(1, getNumFaceElements());
116 case Points:
117 return pair<int,int>(1, 0); //FIXME: dirac
118 default:
119 break;
120 }
121
122 stringstream msg;
123 msg << "getDataShape: Invalid function space type " << fsType
124 << " for " << getDescription();
125 throw RipleyException(msg.str());
126 }
127
128 string RipleyDomain::showTagNames() const
129 {
130 stringstream ret;
131 TagMap::const_iterator it;
132 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
133 if (it!=m_tagMap.begin()) ret << ", ";
134 ret << it->first;
135 }
136 return ret.str();
137 }
138
139 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
140 {
141 /*
142 The idea is to use equivalence classes (i.e. types which can be
143 interpolated back and forth):
144 class 0: DOF <-> Nodes
145 class 1: ReducedDOF <-> ReducedNodes
146 class 2: Points
147 class 3: Elements
148 class 4: ReducedElements
149 class 5: FaceElements
150 class 6: ReducedFaceElements
151
152 There is also a set of lines. Interpolation is possible down a line but not
153 between lines.
154 class 0 and 1 belong to all lines so aren't considered.
155 line 0: class 2
156 line 1: class 3,4
157 line 2: class 5,6
158
159 For classes with multiple members (eg class 1) we have vars to record if
160 there is at least one instance. e.g. hasnodes is true if we have at least
161 one instance of Nodes.
162 */
163 if (fs.empty())
164 return false;
165 vector<bool> hasclass(7, false);
166 vector<int> hasline(3, 0);
167 bool hasnodes=false;
168 bool hasrednodes=false;
169 for (size_t i=0; i<fs.size(); ++i) {
170 switch (fs[i]) {
171 case Nodes: hasnodes=true; // fall through
172 case DegreesOfFreedom:
173 hasclass[0]=true;
174 break;
175 case ReducedNodes: hasrednodes=true; // fall through
176 case ReducedDegreesOfFreedom:
177 hasclass[1]=true;
178 break;
179 case Points:
180 hasline[0]=1;
181 hasclass[2]=true;
182 break;
183 case Elements:
184 hasclass[3]=true;
185 hasline[1]=1;
186 break;
187 case ReducedElements:
188 hasclass[4]=true;
189 hasline[1]=1;
190 break;
191 case FaceElements:
192 hasclass[5]=true;
193 hasline[2]=1;
194 break;
195 case ReducedFaceElements:
196 hasclass[6]=true;
197 hasline[2]=1;
198 break;
199 default:
200 return false;
201 }
202 }
203 int numLines=hasline[0]+hasline[1]+hasline[2];
204
205 // fail if we have more than one leaf group
206 // = there are at least two branches we can't interpolate between
207 if (numLines > 1)
208 return false;
209 else if (numLines==1) {
210 // we have points
211 if (hasline[0]==1)
212 resultcode=Points;
213 else if (hasline[1]==1) {
214 if (hasclass[4])
215 resultcode=ReducedElements;
216 else
217 resultcode=Elements;
218 } else { // hasline[2]==1
219 if (hasclass[6])
220 resultcode=ReducedFaceElements;
221 else
222 resultcode=FaceElements;
223 }
224 } else { // numLines==0
225 if (hasclass[1])
226 // something from class 1
227 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228 else // something from class 0
229 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230 }
231 return true;
232 }
233
234 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
235 int fsType_target) const
236 {
237 if (!isValidFunctionSpaceType(fsType_target)) {
238 stringstream msg;
239 msg << "probeInterpolationOnDomain: Invalid function space type "
240 << fsType_target << " for " << getDescription();
241 throw RipleyException(msg.str());
242 }
243
244 switch (fsType_source) {
245 case Nodes:
246 case DegreesOfFreedom:
247 return true;
248 case ReducedNodes:
249 case ReducedDegreesOfFreedom:
250 return (fsType_target != Nodes &&
251 fsType_target != DegreesOfFreedom);
252 case Elements:
253 return (fsType_target==Elements ||
254 fsType_target==ReducedElements);
255 case FaceElements:
256 return (fsType_target==FaceElements ||
257 fsType_target==ReducedFaceElements);
258 case ReducedElements:
259 case ReducedFaceElements:
260 case Points:
261 return (fsType_target==fsType_source);
262
263 default: {
264 stringstream msg;
265 msg << "probeInterpolationOnDomain: Invalid function space type "
266 << fsType_source << " for " << getDescription();
267 throw RipleyException(msg.str());
268 }
269 }
270 }
271
272 void RipleyDomain::interpolateOnDomain(escript::Data& target,
273 const escript::Data& in) const
274 {
275 const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
276 const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
277 if (inDomain != *this)
278 throw RipleyException("Illegal domain of interpolant");
279 if (targetDomain != *this)
280 throw RipleyException("Illegal domain of interpolation target");
281
282 stringstream msg;
283 msg << "interpolateOnDomain() not implemented for function space "
284 << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
285 << " -> "
286 << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
287
288 const int inFS = in.getFunctionSpace().getTypeCode();
289 const int outFS = target.getFunctionSpace().getTypeCode();
290
291 // simplest case: 1:1 copy
292 if (inFS==outFS) {
293 copyData(target, *const_cast<escript::Data*>(&in));
294 // not allowed: reduced->non-reduced
295 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
296 && (outFS==Nodes || outFS==DegreesOfFreedom)) {
297 throw RipleyException("interpolateOnDomain: Cannot interpolate reduced data to non-reduced data.");
298 } else if ((inFS==Elements && outFS==ReducedElements)
299 || (inFS==FaceElements && outFS==ReducedFaceElements)) {
300 if (in.actsExpanded())
301 averageData(target, *const_cast<escript::Data*>(&in));
302 else
303 copyData(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 case ReducedNodes: //FIXME: reduced
319 copyData(target, *const_cast<escript::Data*>(&in));
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: invalid 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: invalid 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: invalid function space type " << fsType;
567 throw RipleyException(msg.str());
568 }
569 }
570 if (target->size() != num) {
571 target->assign(num, -1);
572 }
573
574 escript::Data& mask=*const_cast<escript::Data*>(&cMask);
575 #pragma omp parallel for
576 for (index_t i=0; i<num; i++) {
577 if (mask.getSampleDataRO(i)[0] > 0) {
578 (*target)[i]=newTag;
579 }
580 }
581 updateTagsInUse(fsType);
582 }
583
584 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
585 {
586 switch(fsType) {
587 case Nodes:
588 if (m_nodeTags.size() > sampleNo)
589 return m_nodeTags[sampleNo];
590 break;
591 case Elements:
592 case ReducedElements:
593 if (m_elementTags.size() > sampleNo)
594 return m_elementTags[sampleNo];
595 break;
596 case FaceElements:
597 case ReducedFaceElements:
598 if (m_faceTags.size() > sampleNo)
599 return m_faceTags[sampleNo];
600 break;
601 default: {
602 stringstream msg;
603 msg << "getTagFromSampleNo: invalid function space type " << fsType;
604 throw RipleyException(msg.str());
605 }
606 }
607 return -1;
608 }
609
610 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
611 {
612 switch(fsType) {
613 case Nodes:
614 return m_nodeTagsInUse.size();
615 case Elements:
616 case ReducedElements:
617 return m_elementTagsInUse.size();
618 case FaceElements:
619 case ReducedFaceElements:
620 return m_faceTagsInUse.size();
621 default: {
622 stringstream msg;
623 msg << "getNumberOfTagsInUse: invalid function space type "
624 << fsType;
625 throw RipleyException(msg.str());
626 }
627 }
628 }
629
630 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
631 {
632 switch(fsType) {
633 case Nodes:
634 return &m_nodeTagsInUse[0];
635 case Elements:
636 case ReducedElements:
637 return &m_elementTagsInUse[0];
638 case FaceElements:
639 case ReducedFaceElements:
640 return &m_faceTagsInUse[0];
641 default: {
642 stringstream msg;
643 msg << "borrowListOfTagsInUse: invalid function space type "
644 << fsType;
645 throw RipleyException(msg.str());
646 }
647 }
648 }
649
650 void RipleyDomain::Print_Mesh_Info(const bool full) const
651 {
652 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
653 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
654 cout << "Number of dimensions: " << m_numDim << endl;
655
656 // write tags
657 if (m_tagMap.size() > 0) {
658 cout << "Tags:" << endl;
659 TagMap::const_iterator it;
660 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
661 cout << " " << setw(5) << it->second << " "
662 << it->first << endl;
663 }
664 }
665 }
666
667 int RipleyDomain::getSystemMatrixTypeId(const int solver,
668 const int preconditioner, const int package, const bool symmetry) const
669 {
670 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
671 package, symmetry, m_mpiInfo);
672 }
673
674 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
675 const int package, const bool symmetry) const
676 {
677 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
678 package, symmetry, m_mpiInfo);
679 }
680
681 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
682 const escript::FunctionSpace& row_functionspace,
683 const int column_blocksize,
684 const escript::FunctionSpace& column_functionspace,
685 const int type) const
686 {
687 bool reduceRowOrder=false;
688 bool reduceColOrder=false;
689 // is the domain right?
690 const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
691 if (row_domain != *this)
692 throw RipleyException("newSystemMatrix: domain of row function space does not match the domain of matrix generator");
693 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
694 if (col_domain != *this)
695 throw RipleyException("newSystemMatrix: domain of column function space does not match the domain of matrix generator");
696 // is the function space type right?
697 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
698 reduceRowOrder=true;
699 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
700 throw RipleyException("newSystemMatrix: illegal function space type for system matrix rows");
701 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
702 reduceColOrder=true;
703 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
704 throw RipleyException("newSystemMatrix: illegal function space type for system matrix columns");
705
706 // generate matrix
707 Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
708 Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
709 row_blocksize, column_blocksize, FALSE);
710 paso::checkPasoError();
711 escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
712 row_functionspace, column_blocksize, column_functionspace));
713 return sma;
714 }
715
716 void RipleyDomain::addPDEToSystem(
717 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
718 const escript::Data& A, const escript::Data& B, const escript::Data& C,
719 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
720 const escript::Data& d, const escript::Data& y,
721 const escript::Data& d_contact, const escript::Data& y_contact,
722 const escript::Data& d_dirac,const escript::Data& y_dirac) const
723 {
724 if (!d_contact.isEmpty() || !y_contact.isEmpty())
725 throw RipleyException("addPDEToSystem: Ripley does not support contact elements");
726
727 paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
728 if (!sma)
729 throw RipleyException("addPDEToSystem: Ripley only accepts Paso system matrices");
730
731 if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
732 throw RipleyException("addPDEToSystem: right hand side coefficients are provided but no right hand side vector given");
733
734 int fsType=UNKNOWN;
735 fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
736 fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
737 fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
738 fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
739 fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
740 fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
741 fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
742 fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
743 fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
744 fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
745 if (fsType==UNKNOWN)
746 return;
747
748 const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
749
750 //TODO: more input param checks (shape, etc)
751
752 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
753
754 if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
755 throw RipleyException("addPDEToSystem: matrix row block size and number of components of right hand side don't match");
756
757 const int numEq=S->logical_row_block_size;
758 const int numComp=S->logical_col_block_size;
759
760 if (numEq != numComp)
761 throw RipleyException("addPDEToSystem: number of equations and number of solutions don't match");
762 //TODO: more system matrix checks
763
764 if (numEq==1) {
765 if (reducedOrder) {
766 assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y);
767 assemblePDEBoundarySingleReduced(S, rhs, d, y);
768 } else {
769 assemblePDESingle(S, rhs, A, B, C, D, X, Y);
770 assemblePDEBoundarySingle(S, rhs, d, y);
771 }
772 } else {
773 if (reducedOrder) {
774 assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y);
775 assemblePDEBoundarySystemReduced(S, rhs, d, y);
776 } else {
777 assemblePDESystem(S, rhs, A, B, C, D, X, Y);
778 assemblePDEBoundarySystem(S, rhs, d, y);
779 }
780 }
781 }
782
783 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
784 const escript::Data& Y, const escript::Data& y,
785 const escript::Data& y_contact, const escript::Data& y_dirac) const
786 {
787 if (!y_contact.isEmpty())
788 throw RipleyException("addPDEToRHS: Ripley does not support contact elements");
789
790 if (rhs.isEmpty()) {
791 if (!X.isEmpty() || !Y.isEmpty())
792 throw RipleyException("addPDEToRHS: right hand side coefficients are provided but no right hand side vector given");
793 else
794 return;
795 }
796
797 if (rhs.getDataPointSize() == 1) {
798 assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
799 assemblePDEBoundarySingle(NULL, rhs, escript::Data(), y);
800 } else {
801 assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
802 assemblePDEBoundarySystem(NULL, rhs, escript::Data(), y);
803 }
804 }
805
806 void RipleyDomain::setNewX(const escript::Data& arg)
807 {
808 throw RipleyException("setNewX(): operation not supported");
809 }
810
811 //protected
812 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
813 {
814 const dim_t numComp = in.getDataPointSize();
815 const dim_t dpp = in.getNumDataPointsPerSample();
816 out.requireWrite();
817 #pragma omp parallel for
818 for (index_t i=0; i<in.getNumSamples(); i++) {
819 const double* src = in.getSampleDataRO(i);
820 double* dest = out.getSampleDataRW(i);
821 for (index_t c=0; c<numComp; c++) {
822 double res=0.;
823 for (index_t q=0; q<dpp; q++)
824 res+=src[c+q*numComp];
825 *dest++ = res/dpp;
826 }
827 }
828 }
829
830 //protected
831 void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
832 {
833 const dim_t numComp = in.getDataPointSize();
834 out.requireWrite();
835 #pragma omp parallel for
836 for (index_t i=0; i<in.getNumSamples(); i++) {
837 const double* src = in.getSampleDataRO(i);
838 copy(src, src+numComp, out.getSampleDataRW(i));
839 }
840 }
841
842 //protected
843 void RipleyDomain::updateTagsInUse(int fsType) const
844 {
845 IndexVector* tagsInUse=NULL;
846 const IndexVector* tags=NULL;
847 switch(fsType) {
848 case Nodes:
849 tags=&m_nodeTags;
850 tagsInUse=&m_nodeTagsInUse;
851 break;
852 case Elements:
853 case ReducedElements:
854 tags=&m_elementTags;
855 tagsInUse=&m_elementTagsInUse;
856 break;
857 case FaceElements:
858 case ReducedFaceElements:
859 tags=&m_faceTags;
860 tagsInUse=&m_faceTagsInUse;
861 break;
862 default:
863 return;
864 }
865
866 // gather global unique tag values from tags into tagsInUse
867 tagsInUse->clear();
868 index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
869
870 while (true) {
871 // find smallest value bigger than lastFoundValue
872 minFoundValue = INDEX_T_MAX;
873 #pragma omp parallel private(local_minFoundValue)
874 {
875 local_minFoundValue = minFoundValue;
876 #pragma omp for schedule(static) nowait
877 for (size_t i = 0; i < tags->size(); i++) {
878 const index_t v = (*tags)[i];
879 if ((v > lastFoundValue) && (v < local_minFoundValue))
880 local_minFoundValue = v;
881 }
882 #pragma omp critical
883 {
884 if (local_minFoundValue < minFoundValue)
885 minFoundValue = local_minFoundValue;
886 }
887 }
888 #ifdef ESYS_MPI
889 local_minFoundValue = minFoundValue;
890 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
891 #endif
892
893 // if we found a new value add it to the tagsInUse vector
894 if (minFoundValue < INDEX_T_MAX) {
895 tagsInUse->push_back(minFoundValue);
896 lastFoundValue = minFoundValue;
897 } else
898 break;
899 }
900 }
901
902 //protected
903 Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
904 const IndexVector& index, const dim_t M, const dim_t N) const
905 {
906 // paso will manage the memory
907 index_t* indexC = MEMALLOC(index.size(), index_t);
908 index_t* ptrC = MEMALLOC(ptr.size(), index_t);
909 copy(index.begin(), index.end(), indexC);
910 copy(ptr.begin(), ptr.end(), ptrC);
911 return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
912 }
913
914 //protected
915 Paso_Pattern* RipleyDomain::createMainPattern() const
916 {
917 IndexVector ptr(1,0);
918 IndexVector index;
919
920 for (index_t i=0; i<getNumDOF(); i++) {
921 // add the DOF itself
922 index.push_back(i);
923 const dim_t num=insertNeighbourNodes(index, i);
924 ptr.push_back(ptr.back()+num+1);
925 }
926
927 return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
928 }
929
930 //protected
931 void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
932 const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
933 {
934 IndexVector ptr(1,0);
935 IndexVector index;
936 for (index_t i=0; i<getNumDOF(); i++) {
937 index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
938 ptr.push_back(ptr.back()+colIndices[i].size());
939 }
940
941 const dim_t M=ptr.size()-1;
942 *colPattern=createPasoPattern(ptr, index, M, N);
943
944 IndexVector rowPtr(1,0);
945 IndexVector rowIndex;
946 for (dim_t id=0; id<N; id++) {
947 dim_t n=0;
948 for (dim_t i=0; i<M; i++) {
949 for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
950 if (index[j]==id) {
951 rowIndex.push_back(i);
952 n++;
953 break;
954 }
955 }
956 }
957 rowPtr.push_back(rowPtr.back()+n);
958 }
959
960 // M and N are now reversed
961 *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
962 }
963
964 //protected
965 void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
966 const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
967 dim_t num_Sol, const vector<double>& array) const
968 {
969 const dim_t numMyCols = mat->pattern->mainPattern->numInput;
970 const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
971 const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
972 const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
973
974 const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
975 const index_t* mainBlock_index = mat->mainBlock->pattern->index;
976 double* mainBlock_val = mat->mainBlock->val;
977 const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
978 const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
979 double* col_coupleBlock_val = mat->col_coupleBlock->val;
980 const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
981 const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
982 double* row_coupleBlock_val = mat->row_coupleBlock->val;
983
984 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
985 // down columns of array
986 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
987 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
988 // only look at the matrix rows stored on this processor
989 if (i_row < numMyRows) {
990 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
991 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
992 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
993 if (i_col < numMyCols) {
994 for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
995 if (mainBlock_index[k] == i_col) {
996 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
997 const dim_t i_Sol=ic+mat->col_block_size*l_col;
998 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
999 const dim_t i_Eq=ir+mat->row_block_size*l_row;
1000 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())];
1001 }
1002 }
1003 break;
1004 }
1005 }
1006 } else {
1007 for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1008 if (col_coupleBlock_index[k] == i_col - numMyCols) {
1009 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1010 const dim_t i_Sol=ic+mat->col_block_size*l_col;
1011 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1012 const dim_t i_Eq=ir+mat->row_block_size*l_row;
1013 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())];
1014 }
1015 }
1016 break;
1017 }
1018 }
1019 }
1020 }
1021 }
1022 } else {
1023 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1024 // across rows of array
1025 for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1026 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1027 if (i_col < numMyCols) {
1028 for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1029 k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1030 {
1031 if (row_coupleBlock_index[k] == i_col) {
1032 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1033 const dim_t i_Sol=ic+mat->col_block_size*l_col;
1034 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1035 const dim_t i_Eq=ir+mat->row_block_size*l_row;
1036 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())];
1037 }
1038 }
1039 break;
1040 }
1041 }
1042 }
1043 }
1044 }
1045 }
1046 }
1047 }
1048 }
1049
1050 //
1051 // the methods that follow have to be implemented by the subclasses
1052 //
1053
1054 string RipleyDomain::getDescription() const
1055 {
1056 throw RipleyException("getDescription() not implemented");
1057 }
1058
1059 void RipleyDomain::write(const string& filename) const
1060 {
1061 throw RipleyException("write() not implemented");
1062 }
1063
1064 void RipleyDomain::dump(const string& filename) const
1065 {
1066 throw RipleyException("dump() not implemented");
1067 }
1068
1069 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
1070 {
1071 throw RipleyException("borrowSampleReferenceIDs() not implemented");
1072 }
1073
1074 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1075 {
1076 throw RipleyException("interpolateACross() not implemented");
1077 }
1078
1079 bool RipleyDomain::probeInterpolationACross(int fsType_source,
1080 const escript::AbstractDomain&, int fsType_target) const
1081 {
1082 throw RipleyException("probeInterpolationACross() not implemented");
1083 }
1084
1085 void RipleyDomain::setToNormal(escript::Data& normal) const
1086 {
1087 throw RipleyException("setToNormal() not implemented");
1088 }
1089
1090 void RipleyDomain::setToSize(escript::Data& size) const
1091 {
1092 throw RipleyException("setToSize() not implemented");
1093 }
1094
1095 bool RipleyDomain::ownSample(int fsType, index_t id) const
1096 {
1097 throw RipleyException("ownSample() 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(
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 dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1181 {
1182 throw RipleyException("insertNeighbourNodes() not implemented");
1183 }
1184
1185 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1186 {
1187 throw RipleyException("assembleCoordinates() not implemented");
1188 }
1189
1190 void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1191 {
1192 throw RipleyException("assembleGradient() not implemented");
1193 }
1194
1195 void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1196 {
1197 throw RipleyException("assembleIntegrate() not implemented");
1198 }
1199
1200 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1201 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1202 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1203 {
1204 throw RipleyException("assemblePDESingle() not implemented");
1205 }
1206
1207 void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1208 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1209 {
1210 throw RipleyException("assemblePDEBoundarySingle() 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
1217 {
1218 throw RipleyException("assemblePDESingleReduced() not implemented");
1219 }
1220
1221 void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1222 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1223 {
1224 throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1225 }
1226
1227 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1228 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1229 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1230 {
1231 throw RipleyException("assemblePDESystem() not implemented");
1232 }
1233
1234 void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1235 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1236 {
1237 throw RipleyException("assemblePDEBoundarySystem() not implemented");
1238 }
1239
1240 void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1241 escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1242 const escript::Data& C, const escript::Data& D, const escript::Data& X,
1243 const escript::Data& Y) const
1244 {
1245 throw RipleyException("assemblePDESystemReduced() not implemented");
1246 }
1247
1248 void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1249 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1250 {
1251 throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1252 }
1253
1254 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1255 {
1256 throw RipleyException("interpolateNodesOnElements() not implemented");
1257 }
1258
1259 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1260 {
1261 throw RipleyException("interpolateNodesOnFaces() not implemented");
1262 }
1263
1264 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1265 {
1266 throw RipleyException("nodesToDOF() not implemented");
1267 }
1268
1269 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1270 {
1271 throw RipleyException("dofToNodes() not implemented");
1272 }
1273
1274 } // end of namespace ripley
1275

  ViewVC Help
Powered by ViewVC 1.1.26