/[escript]/branches/split/ripley/src/RipleyDomain.cpp
ViewVC logotype

Contents of /branches/split/ripley/src/RipleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4747 - (show annotations)
Thu Mar 13 22:52:45 2014 UTC (4 years, 6 months ago) by jfenwick
File size: 55063 byte(s)
switching begins
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include <ripley/RipleyDomain.h>
18 #include <escript/DataFactory.h>
19 #include <escript/FunctionSpaceFactory.h>
20 #include <pasowrap/SystemMatrixAdapter.h>
21 #include <pasowrap/TransportProblemAdapter.h>
22
23 #include <iomanip>
24
25 using namespace std;
26 using paso::SystemMatrixAdapter;
27 using paso::TransportProblemAdapter;
28
29 namespace ripley {
30
31 bool isNotEmpty(string target, map<string, escript::Data> mapping) {
32 map<string, escript::Data>::iterator i = mapping.find(target);
33 return i != mapping.end() && !i->second.isEmpty();
34 }
35
36 void tupleListToMap(map<string, escript::Data>& mapping,
37 boost::python::list& list) {
38 using boost::python::tuple;
39 using boost::python::extract;
40 for (int i = 0; i < len(list); i++) {
41 if (!extract<tuple>(list[i]).check())
42 throw RipleyException("Passed in list contains objects"
43 " other than tuples");
44 tuple t = extract<tuple>(list[i]);
45 if (len(t) != 2 || !extract<std::string>(t[0]).check() ||
46 !extract<escript::Data>(t[1]).check())
47 throw RipleyException("The passed in list must contain tuples"
48 " of the form (string, escript::data)");
49 mapping[extract<std::string>(t[0])] = extract<escript::Data>(t[1]);
50 }
51 }
52
53 RipleyDomain::RipleyDomain(dim_t dim, escript::SubWorld_ptr p) :
54 m_numDim(dim),
55 m_status(0)
56 {
57 if (p.get()==0)
58 {
59 m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
60 }
61 else
62 {
63 m_mpiInfo = Esys_MPIInfo_alloc(p->getMPI()->comm);
64 }
65 }
66
67 RipleyDomain::~RipleyDomain()
68 {
69 Esys_MPIInfo_free(m_mpiInfo);
70 }
71
72 bool RipleyDomain::operator==(const AbstractDomain& other) const
73 {
74 const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
75 if (o) {
76 return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
77 && m_elementTags==o->m_elementTags
78 && m_faceTags==o->m_faceTags);
79 }
80 return false;
81 }
82
83 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
84 {
85 switch (fsType) {
86 case DegreesOfFreedom:
87 case ReducedDegreesOfFreedom:
88 case Nodes:
89 case ReducedNodes:
90 case Elements:
91 case ReducedElements:
92 case FaceElements:
93 case ReducedFaceElements:
94 case Points:
95 return true;
96 default:
97 break;
98 }
99 return false;
100 }
101
102 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
103 {
104 switch (fsType) {
105 case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
106 case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
107 case Nodes: return "Ripley_Nodes";
108 case ReducedNodes: return "Ripley_ReducedNodes";
109 case Elements: return "Ripley_Elements";
110 case ReducedElements: return "Ripley_ReducedElements";
111 case FaceElements: return "Ripley_FaceElements";
112 case ReducedFaceElements: return "Ripley_ReducedFaceElements";
113 case Points: return "Ripley_Points";
114 default:
115 break;
116 }
117 return "Invalid function space type code";
118 }
119
120 pair<int,int> RipleyDomain::getDataShape(int fsType) const
121 {
122 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
123 switch (fsType) {
124 case Nodes:
125 case ReducedNodes: //FIXME: reduced
126 return pair<int,int>(1, getNumNodes());
127 case DegreesOfFreedom:
128 case ReducedDegreesOfFreedom: //FIXME: reduced
129 return pair<int,int>(1, getNumDOF());
130 case Elements:
131 return pair<int,int>(ptsPerSample, getNumElements());
132 case FaceElements:
133 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
134 case ReducedElements:
135 return pair<int,int>(1, getNumElements());
136 case ReducedFaceElements:
137 return pair<int,int>(1, getNumFaceElements());
138 case Points:
139 return pair<int,int>(1, m_diracPoints.size());
140 default:
141 break;
142 }
143
144 stringstream msg;
145 msg << "getDataShape: Invalid function space type " << fsType
146 << " for " << getDescription();
147 throw RipleyException(msg.str());
148 }
149
150 string RipleyDomain::showTagNames() const
151 {
152 stringstream ret;
153 TagMap::const_iterator it;
154 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
155 if (it!=m_tagMap.begin()) ret << ", ";
156 ret << it->first;
157 }
158 return ret.str();
159 }
160
161 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
162 {
163 /*
164 The idea is to use equivalence classes (i.e. types which can be
165 interpolated back and forth):
166 class 0: DOF <-> Nodes
167 class 1: ReducedDOF <-> ReducedNodes
168 class 2: Points
169 class 3: Elements
170 class 4: ReducedElements
171 class 5: FaceElements
172 class 6: ReducedFaceElements
173
174 There is also a set of lines. Interpolation is possible down a line but not
175 between lines.
176 class 0 and 1 belong to all lines so aren't considered.
177 line 0: class 2
178 line 1: class 3,4
179 line 2: class 5,6
180
181 For classes with multiple members (eg class 1) we have vars to record if
182 there is at least one instance. e.g. hasnodes is true if we have at least
183 one instance of Nodes.
184 */
185 if (fs.empty())
186 return false;
187 vector<bool> hasclass(7, false);
188 vector<int> hasline(3, 0);
189 bool hasnodes=false;
190 bool hasrednodes=false;
191 for (size_t i=0; i<fs.size(); ++i) {
192 switch (fs[i]) {
193 case Nodes: hasnodes=true; // fall through
194 case DegreesOfFreedom:
195 hasclass[0]=true;
196 break;
197 case ReducedNodes: hasrednodes=true; // fall through
198 case ReducedDegreesOfFreedom:
199 hasclass[1]=true;
200 break;
201 case Points:
202 hasline[0]=1;
203 hasclass[2]=true;
204 break;
205 case Elements:
206 hasclass[3]=true;
207 hasline[1]=1;
208 break;
209 case ReducedElements:
210 hasclass[4]=true;
211 hasline[1]=1;
212 break;
213 case FaceElements:
214 hasclass[5]=true;
215 hasline[2]=1;
216 break;
217 case ReducedFaceElements:
218 hasclass[6]=true;
219 hasline[2]=1;
220 break;
221 default:
222 return false;
223 }
224 }
225 int numLines=hasline[0]+hasline[1]+hasline[2];
226
227 // fail if we have more than one leaf group
228 // = there are at least two branches we can't interpolate between
229 if (numLines > 1)
230 return false;
231 else if (numLines==1) {
232 // we have points
233 if (hasline[0]==1)
234 resultcode=Points;
235 else if (hasline[1]==1) {
236 if (hasclass[4])
237 resultcode=ReducedElements;
238 else
239 resultcode=Elements;
240 } else { // hasline[2]==1
241 if (hasclass[6])
242 resultcode=ReducedFaceElements;
243 else
244 resultcode=FaceElements;
245 }
246 } else { // numLines==0
247 if (hasclass[1])
248 // something from class 1
249 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
250 else // something from class 0
251 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
252 }
253 return true;
254 }
255
256 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
257 int fsType_target) const
258 {
259 if (!isValidFunctionSpaceType(fsType_target)) {
260 stringstream msg;
261 msg << "probeInterpolationOnDomain: Invalid function space type "
262 << fsType_target << " for " << getDescription();
263 throw RipleyException(msg.str());
264 }
265
266 switch (fsType_source) {
267 case Nodes:
268 case DegreesOfFreedom:
269 return true;
270 case ReducedNodes:
271 case ReducedDegreesOfFreedom:
272 return (fsType_target != Nodes &&
273 fsType_target != DegreesOfFreedom);
274 case Elements:
275 case ReducedElements:
276 return (fsType_target==Elements ||
277 fsType_target==ReducedElements);
278 case FaceElements:
279 case ReducedFaceElements:
280 return (fsType_target==FaceElements ||
281 fsType_target==ReducedFaceElements);
282 case Points:
283 return (fsType_target==fsType_source);
284
285 default: {
286 stringstream msg;
287 msg << "probeInterpolationOnDomain: Invalid function space type "
288 << fsType_source << " for " << getDescription();
289 throw RipleyException(msg.str());
290 }
291 }
292 }
293
294 signed char RipleyDomain::preferredInterpolationOnDomain(int fsType_source,
295 int fsType_target) const
296 {
297 if (!isValidFunctionSpaceType(fsType_target)) {
298 stringstream msg;
299 msg << "preferredInterpolationOnDomain: Invalid function space type "
300 << fsType_target << " for " << getDescription();
301 throw RipleyException(msg.str());
302 }
303
304 if (fsType_source==fsType_target) {
305 return 1;
306 }
307 // There is a complication here in that Nodes and DegreesOfFreedom
308 // can be interpolated to anything, so we need to handle the
309 // reverse case for them specially
310
311 if ((fsType_target==Nodes) || (fsType_target==DegreesOfFreedom))
312 {
313 return -1; // reverse interpolation
314 }
315
316 switch (fsType_source) {
317 case Nodes:
318 case DegreesOfFreedom:
319 return 1;
320 case ReducedNodes:
321 case ReducedDegreesOfFreedom:
322 return (fsType_target != Nodes &&
323 fsType_target != DegreesOfFreedom)?-1:0;
324 case Elements:
325 return (fsType_target==ReducedElements)?1:0;
326 case ReducedElements:
327 return (fsType_target==Elements)?-1:0;
328 case FaceElements:
329 return (fsType_target==ReducedFaceElements)?1:0;
330 case ReducedFaceElements:
331 return (fsType_target==FaceElements)?-1:0;
332 case Points:
333 return false; // other case caught by the if above
334
335 default: {
336 stringstream msg;
337 msg << "probeInterpolationOnDomain: Invalid function space type "
338 << fsType_source << " for " << getDescription();
339 throw RipleyException(msg.str());
340 }
341 }
342 }
343
344 void RipleyDomain::interpolateOnDomain(escript::Data& target,
345 const escript::Data& in) const
346 {
347 const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
348 const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
349 if (inDomain != *this)
350 throw RipleyException("Illegal domain of interpolant");
351 if (targetDomain != *this)
352 throw RipleyException("Illegal domain of interpolation target");
353
354 stringstream msg;
355 msg << "interpolateOnDomain() not implemented for function space "
356 << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
357 << " -> "
358 << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
359
360 const int inFS = in.getFunctionSpace().getTypeCode();
361 const int outFS = target.getFunctionSpace().getTypeCode();
362
363 // simplest case: 1:1 copy
364 if (inFS==outFS) {
365 copyData(target, in);
366 // not allowed: reduced nodes/dof->non-reduced nodes/dof
367 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
368 && (outFS==Nodes || outFS==DegreesOfFreedom)) {
369 throw RipleyException("interpolateOnDomain: Cannot interpolate reduced data to non-reduced data.");
370 } else if ((inFS==Elements && outFS==ReducedElements)
371 || (inFS==FaceElements && outFS==ReducedFaceElements)) {
372 if (in.actsExpanded())
373 averageData(target, in);
374 else
375 copyData(target, in);
376 } else if ((inFS==ReducedElements && outFS==Elements)
377 || (inFS==ReducedFaceElements && outFS==FaceElements)) {
378 multiplyData(target, in);
379 } else {
380 switch (inFS) {
381 case Nodes:
382 case ReducedNodes: //FIXME: reduced
383 switch (outFS) {
384 case DegreesOfFreedom:
385 case ReducedDegreesOfFreedom: //FIXME: reduced
386 if (getMPISize()==1)
387 copyData(target, in);
388 else
389 nodesToDOF(target, in);
390 break;
391
392 case Nodes:
393 case ReducedNodes: //FIXME: reduced
394 copyData(target, in);
395 break;
396 case Elements:
397 interpolateNodesOnElements(target, in, false);
398 break;
399
400 case ReducedElements:
401 interpolateNodesOnElements(target, in, true);
402 break;
403
404 case FaceElements:
405 interpolateNodesOnFaces(target, in, false);
406 break;
407
408 case ReducedFaceElements:
409 interpolateNodesOnFaces(target, in, true);
410 break;
411 case Points:
412 {
413 const dim_t numComp = in.getDataPointSize();
414 target.requireWrite();
415 #pragma omp parallel for
416 for (int i = 0; i < m_diracPoints.size(); i++) {
417 const double* src = in.getSampleDataRO(m_diracPoints[i].node);
418 copy(src, src+numComp, target.getSampleDataRW(i));
419 }
420 }
421 break;
422 default:
423 throw RipleyException(msg.str());
424 }
425 break;
426
427 case DegreesOfFreedom:
428 case ReducedDegreesOfFreedom: //FIXME: reduced
429 switch (outFS) {
430 case Nodes:
431 case ReducedNodes: //FIXME: reduced
432 if (getMPISize()==1)
433 copyData(target, in);
434 else
435 dofToNodes(target, in);
436 break;
437
438 case DegreesOfFreedom:
439 case ReducedDegreesOfFreedom: //FIXME: reduced
440 copyData(target, in);
441 break;
442
443 case Elements:
444 case ReducedElements:
445 if (getMPISize()==1) {
446 interpolateNodesOnElements(target, in, outFS==ReducedElements);
447 } else {
448 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
449 escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
450 interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
451 }
452 break;
453 case FaceElements:
454 case ReducedFaceElements:
455 if (getMPISize()==1) {
456 interpolateNodesOnFaces(target, in, outFS==ReducedFaceElements);
457 } else {
458 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
459 escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
460 interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
461 }
462 break;
463
464 default:
465 throw RipleyException(msg.str());
466 }
467 break;
468 default:
469 throw RipleyException(msg.str());
470 }
471 }
472 }
473
474 escript::Data RipleyDomain::getX() const
475 {
476 return escript::continuousFunction(*this).getX();
477 }
478
479 escript::Data RipleyDomain::getNormal() const
480 {
481 return escript::functionOnBoundary(*this).getNormal();
482 }
483
484 escript::Data RipleyDomain::getSize() const
485 {
486 return escript::function(*this).getSize();
487 }
488
489 void RipleyDomain::setToX(escript::Data& arg) const
490 {
491 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
492 *(arg.getFunctionSpace().getDomain()));
493 if (argDomain != *this)
494 throw RipleyException("setToX: Illegal domain of data point locations");
495 if (!arg.isExpanded())
496 throw RipleyException("setToX: Expanded Data object expected");
497
498 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
499 assembleCoordinates(arg);
500 } else {
501 // interpolate the result
502 escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
503 assembleCoordinates(contData);
504 interpolateOnDomain(arg, contData);
505 }
506 }
507
508 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
509 {
510 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
511 *(arg.getFunctionSpace().getDomain()));
512 if (argDomain != *this)
513 throw RipleyException("setToGradient: Illegal domain of gradient argument");
514 const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
515 *(grad.getFunctionSpace().getDomain()));
516 if (gradDomain != *this)
517 throw RipleyException("setToGradient: Illegal domain of gradient");
518
519 switch (grad.getFunctionSpace().getTypeCode()) {
520 case Elements:
521 case ReducedElements:
522 case FaceElements:
523 case ReducedFaceElements:
524 break;
525 default: {
526 stringstream msg;
527 msg << "setToGradient: not supported for "
528 << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
529 throw RipleyException(msg.str());
530 }
531 }
532
533 switch (arg.getFunctionSpace().getTypeCode()) {
534 case DegreesOfFreedom:
535 case ReducedDegreesOfFreedom:
536 case Nodes:
537 case ReducedNodes:
538 break;
539 default: {
540 throw RipleyException("setToGradient: only supported for nodal input data");
541 }
542 }
543
544 if (getMPISize()>1) {
545 if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
546 escript::Data contArg(arg, escript::continuousFunction(*this));
547 assembleGradient(grad, contArg);
548 } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
549 escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
550 assembleGradient(grad, contArg);
551 } else {
552 assembleGradient(grad, arg);
553 }
554 } else {
555 assembleGradient(grad, arg);
556 }
557 }
558
559 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
560 {
561 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
562 *(arg.getFunctionSpace().getDomain()));
563 if (argDomain != *this)
564 throw RipleyException("setToIntegrals: illegal domain of integration kernel");
565
566 switch (arg.getFunctionSpace().getTypeCode()) {
567 case Nodes:
568 case ReducedNodes:
569 case DegreesOfFreedom:
570 case ReducedDegreesOfFreedom:
571 {
572 escript::Data funcArg(arg, escript::function(*this));
573 assembleIntegrate(integrals, funcArg);
574 }
575 break;
576 case Elements:
577 case ReducedElements:
578 case FaceElements:
579 case ReducedFaceElements:
580 assembleIntegrate(integrals, arg);
581 break;
582 default: {
583 stringstream msg;
584 msg << "setToIntegrals: not supported for "
585 << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
586 throw RipleyException(msg.str());
587 }
588 }
589
590 }
591
592 bool RipleyDomain::isCellOriented(int fsType) const
593 {
594 switch(fsType) {
595 case Nodes:
596 case ReducedNodes:
597 case DegreesOfFreedom:
598 case ReducedDegreesOfFreedom:
599 return false;
600 case Elements:
601 case FaceElements:
602 case Points:
603 case ReducedElements:
604 case ReducedFaceElements:
605 return true;
606 default:
607 break;
608 }
609 stringstream msg;
610 msg << "isCellOriented: invalid function space type " << fsType
611 << " on " << getDescription();
612 throw RipleyException(msg.str());
613 }
614
615 bool RipleyDomain::canTag(int fsType) const
616 {
617 switch(fsType) {
618 case Nodes:
619 case Elements:
620 case ReducedElements:
621 case FaceElements:
622 case Points:
623 case ReducedFaceElements:
624 return true;
625 case DegreesOfFreedom:
626 case ReducedDegreesOfFreedom:
627 case ReducedNodes:
628 return false;
629 default:
630 break;
631 }
632 stringstream msg;
633 msg << "canTag: invalid function space type " << fsType << " on "
634 << getDescription();
635 throw RipleyException(msg.str());
636 }
637
638 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const
639 {
640 IndexVector* target=NULL;
641 dim_t num=0;
642
643 switch(fsType) {
644 case Nodes:
645 num=getNumNodes();
646 target=&m_nodeTags;
647 break;
648 case Elements:
649 case ReducedElements:
650 num=getNumElements();
651 target=&m_elementTags;
652 break;
653 case FaceElements:
654 case ReducedFaceElements:
655 num=getNumFaceElements();
656 target=&m_faceTags;
657 break;
658 default: {
659 stringstream msg;
660 msg << "setTags: invalid function space type " << fsType;
661 throw RipleyException(msg.str());
662 }
663 }
664 if (target->size() != num) {
665 target->assign(num, -1);
666 }
667
668 #pragma omp parallel for
669 for (index_t i=0; i<num; i++) {
670 if (mask.getSampleDataRO(i)[0] > 0) {
671 (*target)[i]=newTag;
672 }
673 }
674 updateTagsInUse(fsType);
675 }
676
677 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
678 {
679 switch(fsType) {
680 case Nodes:
681 if (m_nodeTags.size() > sampleNo)
682 return m_nodeTags[sampleNo];
683 break;
684 case Elements:
685 case ReducedElements:
686 if (m_elementTags.size() > sampleNo)
687 return m_elementTags[sampleNo];
688 break;
689 case Points:
690 if (m_diracPoints.size() > sampleNo)
691 return m_diracPoints[sampleNo].tag;
692 break;
693 case FaceElements:
694 case ReducedFaceElements:
695 if (m_faceTags.size() > sampleNo)
696 return m_faceTags[sampleNo];
697 break;
698 default: {
699 stringstream msg;
700 msg << "getTagFromSampleNo: invalid function space type " << fsType;
701 throw RipleyException(msg.str());
702 }
703 }
704 return -1;
705 }
706
707 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
708 {
709 switch(fsType) {
710 case Nodes:
711 return m_nodeTagsInUse.size();
712 case Elements:
713 case ReducedElements:
714 return m_elementTagsInUse.size();
715 case FaceElements:
716 case ReducedFaceElements:
717 return m_faceTagsInUse.size();
718 default: {
719 stringstream msg;
720 msg << "getNumberOfTagsInUse: invalid function space type "
721 << fsType;
722 throw RipleyException(msg.str());
723 }
724 }
725 }
726
727 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
728 {
729 switch(fsType) {
730 case Nodes:
731 return &m_nodeTagsInUse[0];
732 case Elements:
733 case ReducedElements:
734 return &m_elementTagsInUse[0];
735 case FaceElements:
736 case ReducedFaceElements:
737 return &m_faceTagsInUse[0];
738 default: {
739 stringstream msg;
740 msg << "borrowListOfTagsInUse: invalid function space type "
741 << fsType;
742 throw RipleyException(msg.str());
743 }
744 }
745 }
746
747 void RipleyDomain::Print_Mesh_Info(const bool full) const
748 {
749 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
750 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
751 cout << "Number of dimensions: " << m_numDim << endl;
752 cout << "Number of elements per rank: " << getNumElements() << endl;
753 // write tags
754 if (m_tagMap.size() > 0) {
755 cout << "Tags:" << endl;
756 TagMap::const_iterator it;
757 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
758 cout << " " << setw(5) << it->second << " "
759 << it->first << endl;
760 }
761 }
762 }
763
764 int RipleyDomain::getSystemMatrixTypeId(const int solver,
765 const int preconditioner, const int package, const bool symmetry) const
766 {
767 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
768 package, symmetry, m_mpiInfo);
769 }
770
771 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
772 const int package, const bool symmetry) const
773 {
774 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
775 package, symmetry, m_mpiInfo);
776 }
777
778 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
779 const escript::FunctionSpace& row_functionspace,
780 const int column_blocksize,
781 const escript::FunctionSpace& column_functionspace,
782 const int type) const
783 {
784 bool reduceRowOrder=false;
785 bool reduceColOrder=false;
786 // is the domain right?
787 const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
788 if (row_domain != *this)
789 throw RipleyException("newSystemMatrix: domain of row function space does not match the domain of matrix generator");
790 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
791 if (col_domain != *this)
792 throw RipleyException("newSystemMatrix: domain of column function space does not match the domain of matrix generator");
793 // is the function space type right?
794 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
795 reduceRowOrder=true;
796 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
797 throw RipleyException("newSystemMatrix: illegal function space type for system matrix rows");
798 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
799 reduceColOrder=true;
800 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
801 throw RipleyException("newSystemMatrix: illegal function space type for system matrix columns");
802
803 // generate matrix
804 Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
805 Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
806 row_blocksize, column_blocksize, FALSE);
807 paso::checkPasoError();
808 escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
809 row_functionspace, column_blocksize, column_functionspace));
810 return sma;
811 }
812
813 void RipleyDomain::addToSystem(
814 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
815 std::map<std::string, escript::Data> coefs) const
816 {
817 if (isNotEmpty("d_contact", coefs) || isNotEmpty("y_contact", coefs))
818 throw RipleyException(
819 "addToSystem: Ripley does not support contact elements");
820
821 paso::SystemMatrixAdapter* sma =
822 dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
823 if (!sma)
824 throw RipleyException(
825 "addToSystem: Ripley only accepts Paso system matrices");
826
827 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
828 assemblePDE(S, rhs, coefs);
829 assemblePDEBoundary(S, rhs, coefs);
830 assemblePDEDirac(S, rhs, coefs);
831 }
832
833 void RipleyDomain::addToSystemFromPython(escript::AbstractSystemMatrix& mat,
834 escript::Data& rhs, boost::python::list data) const
835 {
836 std::map<std::string, escript::Data> mapping;
837 tupleListToMap(mapping, data);
838 addToSystem(mat, rhs, mapping);
839 }
840
841 void RipleyDomain::addPDEToSystem(
842 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
843 const escript::Data& A, const escript::Data& B, const escript::Data& C,
844 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
845 const escript::Data& d, const escript::Data& y,
846 const escript::Data& d_contact, const escript::Data& y_contact,
847 const escript::Data& d_dirac,const escript::Data& y_dirac) const
848 {
849 if (!d_contact.isEmpty() || !y_contact.isEmpty())
850 throw RipleyException("addPDEToSystem: Ripley does not support contact elements");
851
852 paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
853 if (!sma)
854 throw RipleyException("addPDEToSystem: Ripley only accepts Paso system matrices");
855
856 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
857
858 std::map<std::string, escript::Data> coefs;
859 coefs["A"] = A; coefs["B"] = B; coefs["C"] = C; coefs["D"] = D;
860 coefs["X"] = X; coefs["Y"] = Y;
861 assemblePDE(S, rhs, coefs);
862
863 std::map<std::string, escript::Data> boundary;
864 boundary["d"] = d;
865 boundary["y"] = y;
866 assemblePDEBoundary(S, rhs, boundary);
867
868 map<string, escript::Data> dirac;
869 dirac["d_dirac"] = d_dirac;
870 dirac["y_dirac"] = y_dirac;
871 assemblePDEDirac(S, rhs, dirac);
872 }
873
874 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
875 const escript::Data& Y, const escript::Data& y,
876 const escript::Data& y_contact, const escript::Data& y_dirac) const
877 {
878 if (!y_contact.isEmpty())
879 throw RipleyException("addPDEToRHS: Ripley does not support contact elements");
880
881 if (rhs.isEmpty()) {
882 if (!X.isEmpty() || !Y.isEmpty())
883 throw RipleyException("addPDEToRHS: right hand side coefficients are provided but no right hand side vector given");
884 else
885 return;
886 }
887 {
888 std::map<std::string, escript::Data> coefs;
889 coefs["X"] = X; coefs["Y"] = Y;
890 assemblePDE(NULL, rhs, coefs);
891 }
892 std::map<std::string, escript::Data> coefs;
893 coefs["y"] = y;
894 assemblePDEBoundary(NULL, rhs, coefs);
895 map<string, escript::Data> dirac;
896 dirac["y_dirac"] = y_dirac;
897 assemblePDEDirac(NULL, rhs, dirac);
898 }
899
900 void RipleyDomain::setAssemblerFromPython(std::string type,
901 boost::python::list options) {
902 std::map<std::string, escript::Data> mapping;
903 tupleListToMap(mapping, options);
904 setAssembler(type, mapping);
905 }
906
907 void RipleyDomain::addToRHSFromPython(escript::Data& rhs,
908 boost::python::list data) const
909 {
910 std::map<std::string, escript::Data> mapping;
911 tupleListToMap(mapping, data);
912 addToRHS(rhs, mapping);
913 }
914
915 void RipleyDomain::addToRHS(escript::Data& rhs,
916 std::map<std::string, escript::Data> coefs) const
917 {
918 if (isNotEmpty("y_contact", coefs))
919 throw RipleyException(
920 "addPDEToRHS: Ripley does not support contact elements");
921
922 if (rhs.isEmpty()) {
923 if (isNotEmpty("X", coefs) || isNotEmpty("Y", coefs))
924 throw RipleyException(
925 "addPDEToRHS: right hand side coefficients are provided "
926 "but no right hand side vector given");
927 else
928 return;
929 }
930
931 assemblePDE(NULL, rhs, coefs);
932 assemblePDEBoundary(NULL, rhs, coefs);
933 assemblePDEDirac(NULL, rhs, coefs);
934 }
935
936 escript::ATP_ptr RipleyDomain::newTransportProblem(const int blocksize,
937 const escript::FunctionSpace& functionspace, const int type) const
938 {
939 bool reduceOrder=false;
940 // is the domain right?
941 const RipleyDomain& domain=dynamic_cast<const RipleyDomain&>(*(functionspace.getDomain()));
942 if (domain != *this)
943 throw RipleyException("newTransportProblem: domain of function space does not match the domain of transport problem generator");
944 // is the function space type right?
945 if (functionspace.getTypeCode()==ReducedDegreesOfFreedom)
946 reduceOrder=true;
947 else if (functionspace.getTypeCode()!=DegreesOfFreedom)
948 throw RipleyException("newTransportProblem: illegal function space type for transport problem");
949
950 // generate matrix
951 Paso_SystemMatrixPattern* pattern=getPattern(reduceOrder, reduceOrder);
952 Paso_TransportProblem* tp = Paso_TransportProblem_alloc(pattern, blocksize);
953 paso::checkPasoError();
954 escript::ATP_ptr atp(new TransportProblemAdapter(tp, blocksize, functionspace));
955 return atp;
956 }
957
958 void RipleyDomain::addPDEToTransportProblem(
959 escript::AbstractTransportProblem& tp,
960 escript::Data& source, const escript::Data& M,
961 const escript::Data& A, const escript::Data& B, const escript::Data& C,
962 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
963 const escript::Data& d, const escript::Data& y,
964 const escript::Data& d_contact, const escript::Data& y_contact,
965 const escript::Data& d_dirac, const escript::Data& y_dirac) const
966 {
967 if (!d_contact.isEmpty() || !y_contact.isEmpty())
968 throw RipleyException("addPDEToTransportProblem: Ripley does not support contact elements");
969
970 paso::TransportProblemAdapter* tpa=dynamic_cast<paso::TransportProblemAdapter*>(&tp);
971 if (!tpa)
972 throw RipleyException("addPDEToTransportProblem: Ripley only accepts Paso transport problems");
973
974 Paso_TransportProblem* ptp = tpa->getPaso_TransportProblem();
975 std::map<std::string, escript::Data> coefs;
976 coefs["D"] = M;
977 assemblePDE(ptp->mass_matrix, source, coefs);
978 coefs["A"] = A; coefs["B"] = B; coefs["C"] = C; coefs["D"] = D;
979 coefs["X"] = X; coefs["Y"] = Y; coefs["d"] = d; coefs["y"] = y;
980 coefs["d_dirac"] = d_dirac; coefs["y_dirac"] = y_dirac;
981 assemblePDE(ptp->transport_matrix, source, coefs);
982 assemblePDEBoundary(ptp->transport_matrix, source, coefs);
983 assemblePDEDirac(ptp->transport_matrix, source, coefs);
984 }
985
986 void RipleyDomain::setNewX(const escript::Data& arg)
987 {
988 throw RipleyException("setNewX(): operation not supported");
989 }
990
991 //protected
992 void RipleyDomain::copyData(escript::Data& out, const escript::Data& in) const
993 {
994 const dim_t numComp = in.getDataPointSize();
995 out.requireWrite();
996 #pragma omp parallel for
997 for (index_t i=0; i<in.getNumSamples(); i++) {
998 const double* src = in.getSampleDataRO(i);
999 copy(src, src+numComp, out.getSampleDataRW(i));
1000 }
1001 }
1002
1003 //protected
1004 void RipleyDomain::averageData(escript::Data& out, const escript::Data& in) const
1005 {
1006 const dim_t numComp = in.getDataPointSize();
1007 const dim_t dpp = in.getNumDataPointsPerSample();
1008 out.requireWrite();
1009 #pragma omp parallel for
1010 for (index_t i=0; i<in.getNumSamples(); i++) {
1011 const double* src = in.getSampleDataRO(i);
1012 double* dest = out.getSampleDataRW(i);
1013 for (index_t c=0; c<numComp; c++) {
1014 double res=0.;
1015 for (index_t q=0; q<dpp; q++)
1016 res+=src[c+q*numComp];
1017 *dest++ = res/dpp;
1018 }
1019 }
1020 }
1021
1022 //protected
1023 void RipleyDomain::multiplyData(escript::Data& out, const escript::Data& in) const
1024 {
1025 const dim_t numComp = in.getDataPointSize();
1026 const dim_t dpp = out.getNumDataPointsPerSample();
1027 out.requireWrite();
1028 #pragma omp parallel for
1029 for (index_t i=0; i<in.getNumSamples(); i++) {
1030 const double* src = in.getSampleDataRO(i);
1031 double* dest = out.getSampleDataRW(i);
1032 for (index_t c=0; c<numComp; c++) {
1033 for (index_t q=0; q<dpp; q++)
1034 dest[c+q*numComp] = src[c];
1035 }
1036 }
1037 }
1038
1039 //protected
1040 void RipleyDomain::updateTagsInUse(int fsType) const
1041 {
1042 IndexVector* tagsInUse=NULL;
1043 const IndexVector* tags=NULL;
1044 switch(fsType) {
1045 case Nodes:
1046 tags=&m_nodeTags;
1047 tagsInUse=&m_nodeTagsInUse;
1048 break;
1049 case Elements:
1050 case ReducedElements:
1051 tags=&m_elementTags;
1052 tagsInUse=&m_elementTagsInUse;
1053 break;
1054 case FaceElements:
1055 case ReducedFaceElements:
1056 tags=&m_faceTags;
1057 tagsInUse=&m_faceTagsInUse;
1058 break;
1059 case Points:
1060 throw RipleyException("updateTagsInUse for Ripley dirac points "
1061 "not supported");
1062 default:
1063 return;
1064 }
1065
1066 // gather global unique tag values from tags into tagsInUse
1067 tagsInUse->clear();
1068 index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
1069
1070 while (true) {
1071 // find smallest value bigger than lastFoundValue
1072 minFoundValue = INDEX_T_MAX;
1073 #pragma omp parallel private(local_minFoundValue)
1074 {
1075 local_minFoundValue = minFoundValue;
1076 long i; // should be size_t but omp mutter mutter
1077 #pragma omp for schedule(static) private(i) nowait
1078 for (i = 0; i < tags->size(); i++) {
1079 const index_t v = (*tags)[i];
1080 if ((v > lastFoundValue) && (v < local_minFoundValue))
1081 local_minFoundValue = v;
1082 }
1083 #pragma omp critical
1084 {
1085 if (local_minFoundValue < minFoundValue)
1086 minFoundValue = local_minFoundValue;
1087 }
1088 }
1089 #ifdef ESYS_MPI
1090 local_minFoundValue = minFoundValue;
1091 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
1092 #endif
1093
1094 // if we found a new value add it to the tagsInUse vector
1095 if (minFoundValue < INDEX_T_MAX) {
1096 tagsInUse->push_back(minFoundValue);
1097 lastFoundValue = minFoundValue;
1098 } else
1099 break;
1100 }
1101 }
1102
1103 //protected
1104 Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
1105 const IndexVector& index, const dim_t M, const dim_t N) const
1106 {
1107 // paso will manage the memory
1108 index_t* indexC = new index_t[index.size()];
1109 index_t* ptrC = new index_t[ptr.size()];
1110 copy(index.begin(), index.end(), indexC);
1111 copy(ptr.begin(), ptr.end(), ptrC);
1112 return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
1113 }
1114
1115 //protected
1116 Paso_Pattern* RipleyDomain::createMainPattern() const
1117 {
1118 IndexVector ptr(1,0);
1119 IndexVector index;
1120
1121 for (index_t i=0; i<getNumDOF(); i++) {
1122 // add the DOF itself
1123 index.push_back(i);
1124 const dim_t num=insertNeighbourNodes(index, i);
1125 ptr.push_back(ptr.back()+num+1);
1126 }
1127
1128 return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
1129 }
1130
1131 //protected
1132 void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
1133 const vector<IndexVector>& rowIndices,
1134 const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
1135 {
1136 IndexVector ptr(1,0);
1137 IndexVector index;
1138 for (index_t i=0; i<getNumDOF(); i++) {
1139 index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
1140 ptr.push_back(ptr.back()+colIndices[i].size());
1141 }
1142
1143 const dim_t M=ptr.size()-1;
1144 *colPattern=createPasoPattern(ptr, index, M, N);
1145
1146 IndexVector rowPtr(1,0);
1147 IndexVector rowIndex;
1148 for (index_t i=0; i<rowIndices.size(); i++) {
1149 rowIndex.insert(rowIndex.end(), rowIndices[i].begin(), rowIndices[i].end());
1150 rowPtr.push_back(rowPtr.back()+rowIndices[i].size());
1151 }
1152
1153 // M and N are now reversed
1154 *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
1155 }
1156
1157 //protected
1158 void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
1159 const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
1160 dim_t num_Sol, const vector<double>& array) const
1161 {
1162 if (mat->type & MATRIX_FORMAT_TRILINOS_CRS)
1163 throw RipleyException("addToSystemMatrix: TRILINOS_CRS not supported");
1164
1165 const dim_t numMyCols = mat->pattern->mainPattern->numInput;
1166 const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
1167 const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
1168 const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
1169
1170 const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
1171 const index_t* mainBlock_index = mat->mainBlock->pattern->index;
1172 double* mainBlock_val = mat->mainBlock->val;
1173 const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
1174 const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
1175 double* col_coupleBlock_val = mat->col_coupleBlock->val;
1176 const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
1177 const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
1178 double* row_coupleBlock_val = mat->row_coupleBlock->val;
1179 index_t offset=(mat->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
1180
1181 #define UPDATE_BLOCK(VAL) do {\
1182 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {\
1183 const dim_t i_Sol=ic+mat->col_block_size*l_col;\
1184 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {\
1185 const dim_t i_Eq=ir+mat->row_block_size*l_row;\
1186 VAL[k*mat->block_size+ir+mat->row_block_size*ic]\
1187 += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];\
1188 }\
1189 }\
1190 } while(0)
1191
1192 if (mat->type & MATRIX_FORMAT_CSC) {
1193 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1194 // down columns of array
1195 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
1196 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1197 if (i_col < numMyCols) {
1198 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
1199 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
1200 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row+offset;
1201 if (i_row < numMyRows+offset) {
1202 for (dim_t k = mainBlock_ptr[i_col]-offset; k < mainBlock_ptr[i_col+1]-offset; ++k) {
1203 if (mainBlock_index[k] == i_row) {
1204 UPDATE_BLOCK(mainBlock_val);
1205 break;
1206 }
1207 }
1208 } else {
1209 for (dim_t k = col_coupleBlock_ptr[i_col]-offset; k < col_coupleBlock_ptr[i_col+1]-offset; ++k) {
1210 if (row_coupleBlock_index[k] == i_row - numMyRows) {
1211 UPDATE_BLOCK(row_coupleBlock_val);
1212 break;
1213 }
1214 }
1215 }
1216 }
1217 }
1218 } else {
1219 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
1220 // across rows of array
1221 for (dim_t l_row=0; l_row<numSubblocks_Eq; ++l_row) {
1222 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row+offset;
1223 if (i_row < numMyRows+offset) {
1224 for (dim_t k = col_coupleBlock_ptr[i_col-numMyCols]-offset;
1225 k < col_coupleBlock_ptr[i_col-numMyCols+1]-offset; ++k)
1226 {
1227 if (col_coupleBlock_index[k] == i_row) {
1228 UPDATE_BLOCK(col_coupleBlock_val);
1229 break;
1230 }
1231 }
1232 }
1233 }
1234 }
1235 }
1236 }
1237 }
1238 } else {
1239 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
1240 // down columns of array
1241 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
1242 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
1243 // only look at the matrix rows stored on this processor
1244 if (i_row < numMyRows) {
1245 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1246 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
1247 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col+offset;
1248 if (i_col < numMyCols+offset) {
1249 for (dim_t k = mainBlock_ptr[i_row]-offset; k < mainBlock_ptr[i_row+1]-offset; ++k) {
1250 if (mainBlock_index[k] == i_col) {
1251 UPDATE_BLOCK(mainBlock_val);
1252 break;
1253 }
1254 }
1255 } else {
1256 for (dim_t k = col_coupleBlock_ptr[i_row]-offset; k < col_coupleBlock_ptr[i_row+1]-offset; ++k) {
1257 if (col_coupleBlock_index[k] == i_col-numMyCols) {
1258 UPDATE_BLOCK(col_coupleBlock_val);
1259 break;
1260 }
1261 }
1262 }
1263 }
1264 }
1265 } else {
1266 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1267 // across rows of array
1268 for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1269 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col+offset;
1270 if (i_col < numMyCols+offset) {
1271 for (dim_t k = row_coupleBlock_ptr[i_row-numMyRows]-offset;
1272 k < row_coupleBlock_ptr[i_row-numMyRows+1]-offset; ++k)
1273 {
1274 if (row_coupleBlock_index[k] == i_col) {
1275 UPDATE_BLOCK(row_coupleBlock_val);
1276 break;
1277 }
1278 }
1279 }
1280 }
1281 }
1282 }
1283 }
1284 }
1285 }
1286 #undef UPDATE_BLOCK
1287 }
1288
1289 //private
1290 void RipleyDomain::assemblePDE(Paso_SystemMatrix* mat, escript::Data& rhs,
1291 std::map<std::string, escript::Data> coefs) const
1292 {
1293 if (rhs.isEmpty() && isNotEmpty("X", coefs) && isNotEmpty("Y", coefs))
1294 throw RipleyException("assemblePDE: right hand side coefficients are "
1295 "provided but no right hand side vector given");
1296
1297 vector<int> fsTypes;
1298 if (isNotEmpty("A", coefs))
1299 fsTypes.push_back(coefs["A"].getFunctionSpace().getTypeCode());
1300 if (isNotEmpty("B", coefs))
1301 fsTypes.push_back(coefs["B"].getFunctionSpace().getTypeCode());
1302 if (isNotEmpty("C", coefs))
1303 fsTypes.push_back(coefs["C"].getFunctionSpace().getTypeCode());
1304 if (isNotEmpty("D", coefs))
1305 fsTypes.push_back(coefs["D"].getFunctionSpace().getTypeCode());
1306 if (isNotEmpty("X", coefs)) //may not exist for some custom assemblers
1307 fsTypes.push_back(coefs["X"].getFunctionSpace().getTypeCode());
1308 else if (isNotEmpty("du", coefs)) //used for (some) custom assemblers
1309 fsTypes.push_back(coefs["du"].getFunctionSpace().getTypeCode());
1310 if (isNotEmpty("Y", coefs))
1311 fsTypes.push_back(coefs["Y"].getFunctionSpace().getTypeCode());
1312
1313 if (fsTypes.empty()) {
1314 return;
1315 }
1316
1317 int fs=fsTypes[0];
1318 if (fs != Elements && fs != ReducedElements)
1319 throw RipleyException("assemblePDE: illegal function space type for coefficients");
1320
1321 for (vector<int>::const_iterator it=fsTypes.begin()+1; it!=fsTypes.end(); it++) {
1322 if (*it != fs) {
1323 throw RipleyException("assemblePDE: coefficient function spaces don't match");
1324 }
1325 }
1326
1327 int numEq, numComp;
1328 if (!mat) {
1329 if (rhs.isEmpty()) {
1330 numEq=numComp=1;
1331 } else {
1332 numEq=numComp=rhs.getDataPointSize();
1333 }
1334 } else {
1335 if (!rhs.isEmpty() && rhs.getDataPointSize()!=mat->logical_row_block_size)
1336 throw RipleyException("assemblePDE: matrix row block size and number of components of right hand side don't match");
1337 numEq = mat->logical_row_block_size;
1338 numComp = mat->logical_col_block_size;
1339 }
1340
1341 if (numEq != numComp)
1342 throw RipleyException("assemblePDE: number of equations and number of solutions don't match");
1343
1344 //TODO: check shape and num samples of coeffs
1345
1346 if (numEq==1) {
1347 if (fs==ReducedElements) {
1348 assembler->assemblePDESingleReduced(mat, rhs, coefs);
1349 } else {
1350 assembler->assemblePDESingle(mat, rhs, coefs);
1351 }
1352 } else {
1353 if (fs==ReducedElements) {
1354 assembler->assemblePDESystemReduced(mat, rhs, coefs);
1355 } else {
1356 assembler->assemblePDESystem(mat, rhs, coefs);
1357 }
1358 }
1359 }
1360
1361 //private
1362 void RipleyDomain::assemblePDEBoundary(Paso_SystemMatrix* mat,
1363 escript::Data& rhs, std::map<std::string, escript::Data> coefs) const
1364 {
1365 std::map<std::string, escript::Data>::iterator iy = coefs.find("y"),
1366 id = coefs.find("d");
1367 if (rhs.isEmpty() && isNotEmpty("y", coefs))
1368 throw RipleyException("assemblePDEBoundary: y provided but no right hand side vector given");
1369
1370 int fs=-1;
1371 if (isNotEmpty("d", coefs))
1372 fs=id->second.getFunctionSpace().getTypeCode();
1373 if (isNotEmpty("y", coefs)) {
1374 if (fs == -1)
1375 fs = iy->second.getFunctionSpace().getTypeCode();
1376 else if (fs != iy->second.getFunctionSpace().getTypeCode())
1377 throw RipleyException("assemblePDEBoundary: coefficient function spaces don't match");
1378 }
1379 if (fs==-1) {
1380 return;
1381 }
1382
1383 if (fs != FaceElements && fs != ReducedFaceElements)
1384 throw RipleyException("assemblePDEBoundary: illegal function space type for coefficients");
1385
1386 int numEq, numComp;
1387 if (!mat) {
1388 if (rhs.isEmpty()) {
1389 numEq=numComp=1;
1390 } else {
1391 numEq=numComp=rhs.getDataPointSize();
1392 }
1393 } else {
1394 if (!rhs.isEmpty() && rhs.getDataPointSize()!=mat->logical_row_block_size)
1395 throw RipleyException("assemblePDEBoundary: matrix row block size and number of components of right hand side don't match");
1396 numEq = mat->logical_row_block_size;
1397 numComp = mat->logical_col_block_size;
1398 }
1399
1400 if (numEq != numComp)
1401 throw RipleyException("assemblePDEBoundary: number of equations and number of solutions don't match");
1402
1403 //TODO: check shape and num samples of coeffs
1404
1405 if (numEq==1) {
1406 if (fs==ReducedFaceElements)
1407 assembler->assemblePDEBoundarySingleReduced(mat, rhs, coefs);
1408 else
1409 assembler->assemblePDEBoundarySingle(mat, rhs, coefs);
1410 } else {
1411 if (fs==ReducedFaceElements)
1412 assembler->assemblePDEBoundarySystemReduced(mat, rhs, coefs);
1413 else
1414 assembler->assemblePDEBoundarySystem(mat, rhs, coefs);
1415 }
1416 }
1417
1418 void RipleyDomain::assemblePDEDirac(Paso_SystemMatrix* mat,
1419 escript::Data& rhs, std::map<std::string, escript::Data> coefs) const
1420 {
1421 bool yNotEmpty = isNotEmpty("y_dirac", coefs),
1422 dNotEmpty = isNotEmpty("d_dirac", coefs);
1423 escript::Data d = dNotEmpty ? coefs["d_dirac"] : escript::Data(),
1424 y = yNotEmpty ? coefs["y_dirac"] : escript::Data();
1425 if (!(yNotEmpty || dNotEmpty)) {
1426 return;
1427 }
1428 int nEq, nComp;
1429 if (!mat) {
1430 if (rhs.isEmpty()) {
1431 nEq=nComp=1;
1432 } else {
1433 nEq=nComp=rhs.getDataPointSize();
1434 }
1435 } else {
1436 if (!rhs.isEmpty() && rhs.getDataPointSize()!=mat->logical_row_block_size)
1437 throw RipleyException("assemblePDEDirac: matrix row block size "
1438 "and number of components of right hand side don't match");
1439 nEq = mat->logical_row_block_size;
1440 nComp = mat->logical_col_block_size;
1441 }
1442 for (int i = 0; i < m_diracPoints.size(); i++) { //only for this rank
1443 IndexVector rowIndex;
1444 rowIndex.push_back(getDofOfNode(m_diracPoints[i].node));
1445 if (yNotEmpty) {
1446 const double *EM_F = y.getSampleDataRO(i);
1447 double *F_p = rhs.getSampleDataRW(0);
1448 for (index_t eq = 0; eq < nEq; eq++) {
1449 F_p[INDEX2(eq, rowIndex[0], nEq)] += EM_F[INDEX2(eq,i,nEq)];
1450 }
1451 }
1452 if (dNotEmpty) {
1453 const double *EM_S = d.getSampleDataRO(i);
1454 std::vector<double> contents(EM_S,
1455 EM_S+mat->row_block_size*nEq*nComp*rowIndex.size());
1456 addToSystemMatrix(mat, rowIndex, nEq, rowIndex, nComp, contents);
1457 }
1458 }
1459 }
1460
1461 bool RipleyDomain::probeInterpolationACross(int fsType_source,
1462 const escript::AbstractDomain&, int fsType_target) const
1463 {
1464 //TODO
1465 return false;
1466 }
1467
1468 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1469 {
1470 throw RipleyException("interpolateACross() not supported");
1471 }
1472
1473 // Expecting ("gaussian", radius, sigma)
1474 bool RipleyDomain::supportsFilter(const boost::python::tuple& t) const
1475 {
1476 if (len(t)==0) { // so we can handle unfiltered randoms
1477 return true;
1478 }
1479 if (len(t)!=3) {
1480 return false;
1481 }
1482 boost::python::extract<string> ex(t[0]);
1483 if (!ex.check() || (ex()!="gaussian"))
1484 {
1485 return false;
1486 }
1487 if (! boost::python::extract<unsigned int>(t[1]).check())
1488 {
1489 return false;
1490 }
1491 return boost::python::extract<double>(t[2]).check();
1492 }
1493
1494 void RipleyDomain::addPoints(int numPoints, const double* points_ptr,
1495 const int* tags_ptr)
1496 {
1497 for (int i = 0; i < numPoints; i++) {
1498 int node = findNode(&points_ptr[i * m_numDim]);
1499 if (node >= 0) {
1500 m_diracPointNodeIDs.push_back(borrowSampleReferenceIDs(Nodes)[node]);
1501 DiracPoint dp;
1502 dp.node = node; //local
1503 dp.tag = tags_ptr[i];
1504 m_diracPoints.push_back(dp);
1505 }
1506 }
1507 }
1508
1509 } // end of namespace ripley
1510

  ViewVC Help
Powered by ViewVC 1.1.26