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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5136 - (show annotations)
Tue Sep 9 07:13:55 2014 UTC (4 years, 5 months ago) by caltinay
File size: 54462 byte(s)
ripley now supports paso solvers again and returns an appropriate matrix type
id. Changed the getSystemMatrixTypeId() method to take a full SolverBuddy
instance and made some other simplifications.

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

  ViewVC Help
Powered by ViewVC 1.1.26