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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 7 months ago) by jfenwick
File size: 55044 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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

  ViewVC Help
Powered by ViewVC 1.1.26