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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5117 - (show annotations)
Mon Aug 18 05:42:09 2014 UTC (4 years, 9 months ago) by caltinay
File size: 51901 byte(s)
Further untangled/simplified paso pattern generation in ripley.

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

  ViewVC Help
Powered by ViewVC 1.1.26