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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6878 - (hide annotations)
Thu Aug 15 06:03:17 2019 UTC (5 weeks, 4 days ago) by aellery
File size: 66046 byte(s)
- Prototype of the EscriptToTVTK function in Weipa. 
- Moved the readthedocs yaml files to the tools folder


1 caltinay 3670
2 jfenwick 3981 /*****************************************************************************
3 caltinay 3670 *
4 jfenwick 6651 * Copyright (c) 2003-2018 by The University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 caltinay 3670 *
7     * Primary Business: Queensland, Australia
8 jfenwick 6112 * Licensed under the Apache License, version 2.0
9     * http://www.apache.org/licenses/LICENSE-2.0
10 caltinay 3670 *
11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 * Development 2012-2013 by School of Earth Sciences
13     * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 *
15     *****************************************************************************/
16 caltinay 3670
17     #include <ripley/RipleyDomain.h>
18 caltinay 5929 #include <ripley/domainhelpers.h>
19    
20 caltinay 3670 #include <escript/DataFactory.h>
21 caltinay 3691 #include <escript/FunctionSpaceFactory.h>
22 caltinay 6169 #include <escript/index.h>
23 caltinay 5929 #include <escript/SolverOptions.h>
24 caltinay 5159
25 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
26 caltinay 6119 #include <trilinoswrap/TrilinosMatrixAdapter.h>
27     #endif
28 caltinay 3670
29 caltinay 6166 #ifdef ESYS_HAVE_PASO
30 caltinay 5929 #include <paso/SystemMatrix.h>
31     #include <paso/Transport.h>
32 caltinay 6166 #endif
33 caltinay 5929
34 caltinay 3670 #include <iomanip>
35 jfenwick 6076 #include <iostream>
36 caltinay 3670
37 caltinay 5120 namespace bp = boost::python;
38    
39 caltinay 3670 using namespace std;
40 caltinay 5984 using escript::ValueError;
41     using escript::NotImplementedError;
42 caltinay 3670
43 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
44 caltinay 6119 using esys_trilinos::TrilinosMatrixAdapter;
45     using esys_trilinos::const_TrilinosGraph_ptr;
46     #endif
47    
48 caltinay 3670 namespace ripley {
49    
50 caltinay 6394 DecompositionPolicy RipleyDomain::m_decompPolicy = DECOMP_ADD_ELEMENTS;
51    
52     void RipleyDomain::setDecompositionPolicy(DecompositionPolicy value)
53     {
54     m_decompPolicy = value;
55     }
56    
57     DecompositionPolicy RipleyDomain::getDecompositionPolicy()
58     {
59     return m_decompPolicy;
60     }
61    
62 caltinay 5120 void tupleListToMap(DataMap& mapping, const bp::list& list)
63 caltinay 5114 {
64 sshaw 4622 for (int i = 0; i < len(list); i++) {
65 caltinay 5120 if (!bp::extract<bp::tuple>(list[i]).check())
66 caltinay 5984 throw ValueError("Passed in list contains objects"
67     " other than tuples");
68 caltinay 5120 bp::tuple t = bp::extract<bp::tuple>(list[i]);
69     if (len(t) != 2 || !bp::extract<string>(t[0]).check() ||
70     !bp::extract<escript::Data>(t[1]).check())
71 caltinay 5984 throw ValueError("The passed in list must contain tuples"
72 caltinay 5148 " of the form (string, escript::Data)");
73 caltinay 5120 mapping[bp::extract<string>(t[0])] = bp::extract<escript::Data>(t[1]);
74 sshaw 4622 }
75     }
76    
77 jfenwick 4934 RipleyDomain::RipleyDomain(dim_t dim, escript::SubWorld_ptr p) :
78 caltinay 3691 m_numDim(dim),
79     m_status(0)
80 caltinay 3670 {
81 caltinay 5114 if (p.get() == NULL)
82 caltinay 5997 m_mpiInfo = escript::makeInfo(MPI_COMM_WORLD);
83 jfenwick 4934 else
84 caltinay 5114 m_mpiInfo = p->getMPI();
85    
86 sshaw 4753 assembler_type = DEFAULT_ASSEMBLER;
87 caltinay 3670 }
88    
89 caltinay 3691 RipleyDomain::~RipleyDomain()
90 caltinay 3670 {
91 jfenwick 4934 // cleanup of MPI is dealt with by shared_ptr
92 caltinay 3670 }
93    
94 caltinay 3744 bool RipleyDomain::operator==(const AbstractDomain& other) const
95     {
96     const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
97     if (o) {
98     return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
99     && m_elementTags==o->m_elementTags
100     && m_faceTags==o->m_faceTags);
101     }
102     return false;
103     }
104    
105 caltinay 3691 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
106 caltinay 3670 {
107 caltinay 3691 switch (fsType) {
108 caltinay 3750 case DegreesOfFreedom:
109     case ReducedDegreesOfFreedom:
110 caltinay 3670 case Nodes:
111     case ReducedNodes:
112     case Elements:
113     case ReducedElements:
114     case FaceElements:
115     case ReducedFaceElements:
116     case Points:
117 caltinay 3691 return true;
118     default:
119 caltinay 3670 break;
120     }
121 caltinay 3691 return false;
122 caltinay 3670 }
123    
124 caltinay 3691 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
125 caltinay 3670 {
126 caltinay 3691 switch (fsType) {
127 caltinay 5120 case DegreesOfFreedom:
128     return "Ripley_DegreesOfFreedom [Solution(domain)]";
129     case ReducedDegreesOfFreedom:
130     return "Ripley_ReducedDegreesOfFreedom [ReducedSolution(domain)]";
131     case Nodes:
132     return "Ripley_Nodes [ContinuousFunction(domain)]";
133     case ReducedNodes:
134     return "Ripley_ReducedNodes [ReducedContinuousFunction(domain)]";
135     case Elements:
136     return "Ripley_Elements [Function(domain)]";
137     case ReducedElements:
138     return "Ripley_ReducedElements [ReducedFunction(domain)]";
139     case FaceElements:
140     return "Ripley_FaceElements [FunctionOnBoundary(domain)]";
141     case ReducedFaceElements:
142     return "Ripley_ReducedFaceElements [ReducedFunctionOnBoundary(domain)]";
143     case Points:
144     return "Ripley_Points [DiracDeltaFunctions(domain)]";
145 caltinay 3691 default:
146 caltinay 3670 break;
147     }
148 caltinay 3691 return "Invalid function space type code";
149 caltinay 3670 }
150    
151 caltinay 5120 pair<int,dim_t> RipleyDomain::getDataShape(int fsType) const
152 caltinay 3670 {
153 caltinay 3691 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
154     switch (fsType) {
155 caltinay 3670 case Nodes:
156 caltinay 5120 case ReducedNodes:
157     return pair<int,dim_t>(1, getNumNodes());
158 caltinay 3750 case DegreesOfFreedom:
159 caltinay 5120 case ReducedDegreesOfFreedom:
160     return pair<int,dim_t>(1, getNumDOF());
161 caltinay 3670 case Elements:
162 caltinay 5120 return pair<int,dim_t>(ptsPerSample, getNumElements());
163 caltinay 3670 case FaceElements:
164 caltinay 5120 return pair<int,dim_t>(ptsPerSample/2, getNumFaceElements());
165 caltinay 3711 case ReducedElements:
166 caltinay 5120 return pair<int,dim_t>(1, getNumElements());
167 caltinay 3711 case ReducedFaceElements:
168 caltinay 5120 return pair<int,dim_t>(1, getNumFaceElements());
169 caltinay 3697 case Points:
170 caltinay 5120 return pair<int,dim_t>(1, m_diracPoints.size());
171 caltinay 3670 default:
172     break;
173     }
174    
175 caltinay 3691 stringstream msg;
176 caltinay 3791 msg << "getDataShape: Invalid function space type " << fsType
177 caltinay 3740 << " for " << getDescription();
178 caltinay 5984 throw ValueError(msg.str());
179 caltinay 3670 }
180    
181 caltinay 3691 string RipleyDomain::showTagNames() const
182 caltinay 3670 {
183 caltinay 3691 stringstream ret;
184     TagMap::const_iterator it;
185     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
186     if (it!=m_tagMap.begin()) ret << ", ";
187     ret << it->first;
188 caltinay 3670 }
189 caltinay 3691 return ret.str();
190 caltinay 3670 }
191    
192     bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
193     {
194 caltinay 3691 /*
195     The idea is to use equivalence classes (i.e. types which can be
196     interpolated back and forth):
197 caltinay 3750 class 0: DOF <-> Nodes
198     class 1: ReducedDOF <-> ReducedNodes
199 caltinay 3740 class 2: Points
200     class 3: Elements
201     class 4: ReducedElements
202     class 5: FaceElements
203     class 6: ReducedFaceElements
204 caltinay 3670
205 caltinay 3691 There is also a set of lines. Interpolation is possible down a line but not
206     between lines.
207 caltinay 3740 class 0 and 1 belong to all lines so aren't considered.
208     line 0: class 2
209     line 1: class 3,4
210     line 2: class 5,6
211 caltinay 3750
212     For classes with multiple members (eg class 1) we have vars to record if
213     there is at least one instance. e.g. hasnodes is true if we have at least
214     one instance of Nodes.
215 caltinay 3670 */
216     if (fs.empty())
217     return false;
218 caltinay 3740 vector<bool> hasclass(7, false);
219     vector<int> hasline(3, 0);
220 caltinay 3750 bool hasnodes=false;
221     bool hasrednodes=false;
222 caltinay 3740 for (size_t i=0; i<fs.size(); ++i) {
223 caltinay 3670 switch (fs[i]) {
224 caltinay 3750 case Nodes: hasnodes=true; // fall through
225     case DegreesOfFreedom:
226 caltinay 3740 hasclass[0]=true;
227 caltinay 3670 break;
228 caltinay 3750 case ReducedNodes: hasrednodes=true; // fall through
229     case ReducedDegreesOfFreedom:
230 caltinay 3740 hasclass[1]=true;
231 caltinay 3670 break;
232     case Points:
233     hasline[0]=1;
234 caltinay 3740 hasclass[2]=true;
235 caltinay 3670 break;
236     case Elements:
237 caltinay 3750 hasclass[3]=true;
238 caltinay 3670 hasline[1]=1;
239     break;
240     case ReducedElements:
241 caltinay 3740 hasclass[4]=true;
242 caltinay 3670 hasline[1]=1;
243     break;
244     case FaceElements:
245 caltinay 3750 hasclass[5]=true;
246 caltinay 3670 hasline[2]=1;
247     break;
248     case ReducedFaceElements:
249 caltinay 3740 hasclass[6]=true;
250 caltinay 3670 hasline[2]=1;
251     break;
252     default:
253     return false;
254     }
255     }
256 caltinay 3740 int numLines=hasline[0]+hasline[1]+hasline[2];
257 caltinay 3670
258     // fail if we have more than one leaf group
259     // = there are at least two branches we can't interpolate between
260 caltinay 3740 if (numLines > 1)
261 caltinay 3670 return false;
262 caltinay 3740 else if (numLines==1) {
263 caltinay 3750 // we have points
264 caltinay 3691 if (hasline[0]==1)
265 caltinay 3670 resultcode=Points;
266 caltinay 3691 else if (hasline[1]==1) {
267 caltinay 3740 if (hasclass[4])
268 caltinay 3670 resultcode=ReducedElements;
269 caltinay 3691 else
270 caltinay 3670 resultcode=Elements;
271 caltinay 3740 } else { // hasline[2]==1
272     if (hasclass[6])
273 caltinay 3670 resultcode=ReducedFaceElements;
274 caltinay 3691 else
275 caltinay 3670 resultcode=FaceElements;
276 caltinay 3740 }
277     } else { // numLines==0
278     if (hasclass[1])
279 caltinay 3750 // something from class 1
280     resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
281     else // something from class 0
282     resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
283 caltinay 3670 }
284     return true;
285     }
286    
287 caltinay 3732 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
288     int fsType_target) const
289     {
290 caltinay 3764 if (!isValidFunctionSpaceType(fsType_target)) {
291 caltinay 3732 stringstream msg;
292 caltinay 3791 msg << "probeInterpolationOnDomain: Invalid function space type "
293 caltinay 3740 << fsType_target << " for " << getDescription();
294 caltinay 5984 throw ValueError(msg.str());
295 caltinay 3732 }
296    
297     switch (fsType_source) {
298     case Nodes:
299 caltinay 3750 case DegreesOfFreedom:
300 caltinay 3732 return true;
301     case ReducedNodes:
302 caltinay 3750 case ReducedDegreesOfFreedom:
303     return (fsType_target != Nodes &&
304     fsType_target != DegreesOfFreedom);
305 caltinay 3732 case Elements:
306 caltinay 4211 case ReducedElements:
307 caltinay 3732 return (fsType_target==Elements ||
308     fsType_target==ReducedElements);
309     case FaceElements:
310 caltinay 4211 case ReducedFaceElements:
311 caltinay 3732 return (fsType_target==FaceElements ||
312     fsType_target==ReducedFaceElements);
313     case Points:
314 caltinay 3764 return (fsType_target==fsType_source);
315 caltinay 3732
316     default: {
317     stringstream msg;
318 caltinay 3791 msg << "probeInterpolationOnDomain: Invalid function space type "
319 caltinay 3740 << fsType_source << " for " << getDescription();
320 caltinay 5984 throw ValueError(msg.str());
321 caltinay 3732 }
322     }
323     }
324    
325 jfenwick 4255 signed char RipleyDomain::preferredInterpolationOnDomain(int fsType_source,
326 caltinay 5120 int fsType_target) const
327 jfenwick 4255 {
328     if (!isValidFunctionSpaceType(fsType_target)) {
329     stringstream msg;
330     msg << "preferredInterpolationOnDomain: Invalid function space type "
331     << fsType_target << " for " << getDescription();
332 caltinay 5984 throw ValueError(msg.str());
333 jfenwick 4255 }
334    
335     if (fsType_source==fsType_target) {
336     return 1;
337     }
338     // There is a complication here in that Nodes and DegreesOfFreedom
339     // can be interpolated to anything, so we need to handle the
340     // reverse case for them specially
341    
342 caltinay 5120 if ((fsType_target==Nodes) || (fsType_target==DegreesOfFreedom)) {
343     return -1; // reverse interpolation
344 jfenwick 4255 }
345    
346     switch (fsType_source) {
347     case Nodes:
348     case DegreesOfFreedom:
349     return 1;
350     case ReducedNodes:
351     case ReducedDegreesOfFreedom:
352     return (fsType_target != Nodes &&
353 caltinay 5120 fsType_target != DegreesOfFreedom) ? -1 : 0;
354 jfenwick 4255 case Elements:
355 caltinay 5120 return (fsType_target==ReducedElements) ? 1 : 0;
356 jfenwick 4255 case ReducedElements:
357 caltinay 5120 return (fsType_target==Elements) ? -1 : 0;
358 jfenwick 4255 case FaceElements:
359 caltinay 5120 return (fsType_target==ReducedFaceElements) ? 1 : 0;
360 jfenwick 4255 case ReducedFaceElements:
361 caltinay 5120 return (fsType_target==FaceElements) ? -1 : 0;
362     case Points:
363     return false; // other case caught by the if above
364 jfenwick 4255 default: {
365     stringstream msg;
366     msg << "probeInterpolationOnDomain: Invalid function space type "
367     << fsType_source << " for " << getDescription();
368 caltinay 5984 throw ValueError(msg.str());
369 jfenwick 4255 }
370     }
371     }
372    
373 caltinay 3701 void RipleyDomain::interpolateOnDomain(escript::Data& target,
374     const escript::Data& in) const
375     {
376     const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
377     const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
378     if (inDomain != *this)
379 caltinay 5984 throw ValueError("Illegal domain of interpolant");
380 caltinay 3701 if (targetDomain != *this)
381 caltinay 5984 throw ValueError("Illegal domain of interpolation target");
382 caltinay 6435 if (target.isComplex() != in.isComplex())
383     throw ValueError("Complexity of input and output must match");
384 caltinay 3701
385     stringstream msg;
386     msg << "interpolateOnDomain() not implemented for function space "
387     << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
388     << " -> "
389     << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
390    
391 caltinay 3744 const int inFS = in.getFunctionSpace().getTypeCode();
392     const int outFS = target.getFunctionSpace().getTypeCode();
393 caltinay 3701
394 caltinay 3744 // simplest case: 1:1 copy
395     if (inFS==outFS) {
396 caltinay 6435 if (in.isComplex())
397     copyData<cplx_t>(target, in);
398     else
399     copyData<real_t>(target, in);
400 caltinay 5120 // not allowed: reduced nodes/DOF->non-reduced nodes/DOF
401 caltinay 3750 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
402     && (outFS==Nodes || outFS==DegreesOfFreedom)) {
403 caltinay 5984 throw ValueError("interpolateOnDomain: Cannot interpolate "
404 caltinay 5120 "reduced data to non-reduced data.");
405 caltinay 3744 } else if ((inFS==Elements && outFS==ReducedElements)
406     || (inFS==FaceElements && outFS==ReducedFaceElements)) {
407 caltinay 6435 if (in.actsExpanded()) {
408     if (in.isComplex())
409     averageData<cplx_t>(target, in);
410     else
411     averageData<real_t>(target, in);
412     } else {
413     if (in.isComplex())
414     copyData<cplx_t>(target, in);
415     else
416     copyData<real_t>(target, in);
417     }
418 caltinay 4211 } else if ((inFS==ReducedElements && outFS==Elements)
419     || (inFS==ReducedFaceElements && outFS==FaceElements)) {
420 caltinay 6435 if (in.isComplex())
421     multiplyData<cplx_t>(target, in);
422     else
423     multiplyData<real_t>(target, in);
424 caltinay 3744 } else {
425     switch (inFS) {
426     case Nodes:
427 caltinay 5120 case ReducedNodes:
428 caltinay 3744 switch (outFS) {
429 caltinay 3750 case DegreesOfFreedom:
430 caltinay 5120 case ReducedDegreesOfFreedom:
431 caltinay 6580 if (getMPISize()==1) {
432 caltinay 6435 if (in.isComplex())
433     copyData<cplx_t>(target, in);
434     else
435     copyData<real_t>(target, in);
436 caltinay 6580 } else {
437 caltinay 6581 if (in.isComplex())
438     throw NotImplementedError("nodesToDOF not implemented for complex Data");
439     else
440     nodesToDOF(target, in);
441 caltinay 6580 }
442 caltinay 3744 break;
443 caltinay 3701
444 caltinay 3755 case Nodes:
445 caltinay 5120 case ReducedNodes:
446 caltinay 6435 if (in.isComplex())
447     copyData<cplx_t>(target, in);
448     else
449     copyData<real_t>(target, in);
450 caltinay 3755 break;
451 caltinay 3744 case Elements:
452 caltinay 4626 interpolateNodesOnElements(target, in, false);
453 caltinay 3744 break;
454 caltinay 3711
455 caltinay 3744 case ReducedElements:
456 caltinay 4626 interpolateNodesOnElements(target, in, true);
457 caltinay 3744 break;
458 caltinay 3701
459 caltinay 3744 case FaceElements:
460 caltinay 4626 interpolateNodesOnFaces(target, in, false);
461 caltinay 3744 break;
462 caltinay 3711
463 caltinay 3744 case ReducedFaceElements:
464 caltinay 4626 interpolateNodesOnFaces(target, in, true);
465 caltinay 3744 break;
466 sshaw 4629 case Points:
467     {
468     const dim_t numComp = in.getDataPointSize();
469 caltinay 5039 const int nDirac = m_diracPoints.size();
470 sshaw 4629 target.requireWrite();
471 caltinay 5120 #pragma omp parallel for
472 caltinay 5039 for (int i = 0; i < nDirac; i++) {
473 sshaw 4629 const double* src = in.getSampleDataRO(m_diracPoints[i].node);
474     copy(src, src+numComp, target.getSampleDataRW(i));
475     }
476     }
477     break;
478 caltinay 3744 default:
479 caltinay 5984 throw NotImplementedError(msg.str());
480 caltinay 3744 }
481     break;
482 caltinay 3750
483     case DegreesOfFreedom:
484 caltinay 5120 case ReducedDegreesOfFreedom:
485 caltinay 3755 switch (outFS) {
486     case Nodes:
487 caltinay 5120 case ReducedNodes:
488 caltinay 3755 if (getMPISize()==1)
489 caltinay 6435 if (in.isComplex())
490     copyData<cplx_t>(target, in);
491     else
492     copyData<real_t>(target, in);
493 caltinay 3755 else
494 caltinay 6580 if (in.isComplex())
495     dofToNodes<cplx_t>(target, in);
496     else
497     dofToNodes<real_t>(target, in);
498 caltinay 3755 break;
499    
500     case DegreesOfFreedom:
501 caltinay 5120 case ReducedDegreesOfFreedom:
502 caltinay 6435 if (in.isComplex())
503     copyData<cplx_t>(target, in);
504     else
505     copyData<real_t>(target, in);
506 caltinay 3755 break;
507    
508     case Elements:
509     case ReducedElements:
510     if (getMPISize()==1) {
511 caltinay 4626 interpolateNodesOnElements(target, in, outFS==ReducedElements);
512 caltinay 3755 } else {
513 caltinay 3764 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
514     escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
515 caltinay 3755 interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
516     }
517     break;
518     case FaceElements:
519     case ReducedFaceElements:
520     if (getMPISize()==1) {
521 caltinay 4626 interpolateNodesOnFaces(target, in, outFS==ReducedFaceElements);
522 caltinay 3755 } else {
523 caltinay 3764 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
524     escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
525     interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
526 caltinay 3755 }
527     break;
528    
529     default:
530 caltinay 5984 throw NotImplementedError(msg.str());
531 caltinay 3755 }
532 caltinay 3754 break;
533 sshaw 5317 case Points:
534     switch(outFS) {
535     case Nodes:
536     {
537     const dim_t numComp = in.getDataPointSize();
538     const int nDirac = m_diracPoints.size();
539     target.requireWrite();
540     #pragma omp parallel for
541     for (int i = 0; i < nDirac; i++) {
542     const double* src = in.getSampleDataRO(i);
543     copy(src, src+numComp, target.getSampleDataRW(m_diracPoints[i].node));
544     }
545     }
546    
547     }
548     break;
549 caltinay 3744 default:
550 caltinay 5984 throw NotImplementedError(msg.str());
551 caltinay 3744 }
552 caltinay 3701 }
553     }
554    
555 caltinay 3691 escript::Data RipleyDomain::getX() const
556 caltinay 3670 {
557 caltinay 3691 return escript::continuousFunction(*this).getX();
558     }
559    
560 aellery 6872 #ifdef ESYS_HAVE_BOOST_NUMPY
561     boost::python::numpy::ndarray RipleyDomain::getNumpyX() const
562     {
563     return escript::continuousFunction(*this).getNumpyX();
564     }
565 aellery 6878
566     boost::python::numpy::ndarray RipleyDomain::getConnectivityInfo() const
567     {
568     throw RipleyException("This feature is currently not supported by Ripley.");
569     }
570 aellery 6872 #endif
571    
572 caltinay 3691 escript::Data RipleyDomain::getNormal() const
573     {
574     return escript::functionOnBoundary(*this).getNormal();
575     }
576    
577     escript::Data RipleyDomain::getSize() const
578     {
579     return escript::function(*this).getSize();
580     }
581    
582     void RipleyDomain::setToX(escript::Data& arg) const
583     {
584     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
585     *(arg.getFunctionSpace().getDomain()));
586     if (argDomain != *this)
587 caltinay 5984 throw ValueError("setToX: Illegal domain of data point locations");
588 caltinay 3691 if (!arg.isExpanded())
589 caltinay 5984 throw ValueError("setToX: Expanded Data object expected");
590 caltinay 3691
591     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
592     assembleCoordinates(arg);
593     } else {
594     // interpolate the result
595     escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
596     assembleCoordinates(contData);
597     interpolateOnDomain(arg, contData);
598 caltinay 3670 }
599 caltinay 3691 }
600 caltinay 3670
601 caltinay 3764 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
602     {
603     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
604     *(arg.getFunctionSpace().getDomain()));
605     if (argDomain != *this)
606 caltinay 5984 throw ValueError("setToGradient: Illegal domain of gradient argument");
607 caltinay 3764 const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
608     *(grad.getFunctionSpace().getDomain()));
609     if (gradDomain != *this)
610 caltinay 5984 throw ValueError("setToGradient: Illegal domain of gradient");
611 caltinay 3764
612     switch (grad.getFunctionSpace().getTypeCode()) {
613     case Elements:
614     case ReducedElements:
615     case FaceElements:
616     case ReducedFaceElements:
617     break;
618     default: {
619     stringstream msg;
620 caltinay 3791 msg << "setToGradient: not supported for "
621 caltinay 3764 << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
622 caltinay 5984 throw ValueError(msg.str());
623 caltinay 3764 }
624     }
625    
626 caltinay 4058 switch (arg.getFunctionSpace().getTypeCode()) {
627     case DegreesOfFreedom:
628     case ReducedDegreesOfFreedom:
629     case Nodes:
630     case ReducedNodes:
631     break;
632     default: {
633 caltinay 5984 throw ValueError("setToGradient: only supported for nodal input data");
634 caltinay 4058 }
635     }
636    
637 caltinay 5120 if (getMPISize() > 1) {
638 caltinay 3764 if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
639     escript::Data contArg(arg, escript::continuousFunction(*this));
640     assembleGradient(grad, contArg);
641     } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
642     escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
643     assembleGradient(grad, contArg);
644     } else {
645 caltinay 4626 assembleGradient(grad, arg);
646 caltinay 3764 }
647     } else {
648 caltinay 4626 assembleGradient(grad, arg);
649 caltinay 3764 }
650     }
651    
652 caltinay 6473 void RipleyDomain::setToIntegrals(vector<real_t>& integrals, const escript::Data& arg) const
653 caltinay 3764 {
654 caltinay 6473 setToIntegralsWorker<real_t>(integrals, arg);
655     }
656    
657     void RipleyDomain::setToIntegrals(vector<cplx_t>& integrals, const escript::Data& arg) const
658     {
659     setToIntegralsWorker<cplx_t>(integrals, arg);
660     }
661    
662     template<typename Scalar>
663     void RipleyDomain::setToIntegralsWorker(std::vector<Scalar>& integrals,
664     const escript::Data& arg) const
665     {
666     const RipleyDomain& argDomain = dynamic_cast<const RipleyDomain&>(
667 caltinay 3764 *(arg.getFunctionSpace().getDomain()));
668     if (argDomain != *this)
669 caltinay 5984 throw ValueError("setToIntegrals: illegal domain of integration kernel");
670 caltinay 3764
671     switch (arg.getFunctionSpace().getTypeCode()) {
672     case Nodes:
673     case ReducedNodes:
674     case DegreesOfFreedom:
675     case ReducedDegreesOfFreedom:
676     {
677     escript::Data funcArg(arg, escript::function(*this));
678     assembleIntegrate(integrals, funcArg);
679     }
680     break;
681     case Elements:
682     case ReducedElements:
683     case FaceElements:
684     case ReducedFaceElements:
685 caltinay 4626 assembleIntegrate(integrals, arg);
686 caltinay 3764 break;
687     default: {
688     stringstream msg;
689 caltinay 3791 msg << "setToIntegrals: not supported for "
690 caltinay 3764 << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
691 caltinay 5984 throw ValueError(msg.str());
692 caltinay 3764 }
693     }
694     }
695    
696 caltinay 3691 bool RipleyDomain::isCellOriented(int fsType) const
697     {
698     switch(fsType) {
699 caltinay 3670 case Nodes:
700 caltinay 3764 case ReducedNodes:
701 caltinay 3750 case DegreesOfFreedom:
702     case ReducedDegreesOfFreedom:
703 caltinay 3691 return false;
704 caltinay 3670 case Elements:
705 caltinay 3691 case FaceElements:
706     case Points:
707 caltinay 3670 case ReducedElements:
708     case ReducedFaceElements:
709 caltinay 3691 return true;
710 caltinay 3670 default:
711 caltinay 3691 break;
712 caltinay 3670 }
713 caltinay 3691 stringstream msg;
714 caltinay 3791 msg << "isCellOriented: invalid function space type " << fsType
715 caltinay 3740 << " on " << getDescription();
716 caltinay 5984 throw ValueError(msg.str());
717 caltinay 3670 }
718    
719 caltinay 3722 bool RipleyDomain::canTag(int fsType) const
720     {
721     switch(fsType) {
722     case Nodes:
723     case Elements:
724     case ReducedElements:
725     case FaceElements:
726 sshaw 4622 case Points:
727 caltinay 3722 case ReducedFaceElements:
728     return true;
729 caltinay 3750 case DegreesOfFreedom:
730     case ReducedDegreesOfFreedom:
731 caltinay 3764 case ReducedNodes:
732 caltinay 3722 return false;
733     default:
734     break;
735     }
736     stringstream msg;
737 caltinay 3791 msg << "canTag: invalid function space type " << fsType << " on "
738 caltinay 3740 << getDescription();
739 caltinay 5984 throw ValueError(msg.str());
740 caltinay 3722 }
741    
742 caltinay 5187 void RipleyDomain::setTags(int fsType, int newTag, const escript::Data& mask) const
743 caltinay 3722 {
744 caltinay 5187 vector<int>* target=NULL;
745 caltinay 3722 dim_t num=0;
746    
747     switch(fsType) {
748     case Nodes:
749     num=getNumNodes();
750     target=&m_nodeTags;
751     break;
752     case Elements:
753     case ReducedElements:
754     num=getNumElements();
755     target=&m_elementTags;
756     break;
757     case FaceElements:
758     case ReducedFaceElements:
759     num=getNumFaceElements();
760     target=&m_faceTags;
761     break;
762     default: {
763     stringstream msg;
764 caltinay 3791 msg << "setTags: invalid function space type " << fsType;
765 caltinay 5984 throw ValueError(msg.str());
766 caltinay 3722 }
767     }
768     if (target->size() != num) {
769     target->assign(num, -1);
770     }
771    
772     #pragma omp parallel for
773     for (index_t i=0; i<num; i++) {
774     if (mask.getSampleDataRO(i)[0] > 0) {
775     (*target)[i]=newTag;
776     }
777     }
778     updateTagsInUse(fsType);
779     }
780    
781 caltinay 5120 int RipleyDomain::getTagFromSampleNo(int fsType, dim_t sampleNo) const
782 caltinay 3722 {
783     switch(fsType) {
784     case Nodes:
785     if (m_nodeTags.size() > sampleNo)
786     return m_nodeTags[sampleNo];
787     break;
788     case Elements:
789     case ReducedElements:
790     if (m_elementTags.size() > sampleNo)
791     return m_elementTags[sampleNo];
792     break;
793 sshaw 4622 case Points:
794 sshaw 4629 if (m_diracPoints.size() > sampleNo)
795     return m_diracPoints[sampleNo].tag;
796 sshaw 4622 break;
797 caltinay 3722 case FaceElements:
798     case ReducedFaceElements:
799     if (m_faceTags.size() > sampleNo)
800     return m_faceTags[sampleNo];
801     break;
802     default: {
803     stringstream msg;
804 caltinay 3791 msg << "getTagFromSampleNo: invalid function space type " << fsType;
805 caltinay 5984 throw ValueError(msg.str());
806 caltinay 3722 }
807     }
808     return -1;
809     }
810    
811     int RipleyDomain::getNumberOfTagsInUse(int fsType) const
812     {
813     switch(fsType) {
814     case Nodes:
815     return m_nodeTagsInUse.size();
816     case Elements:
817     case ReducedElements:
818     return m_elementTagsInUse.size();
819     case FaceElements:
820     case ReducedFaceElements:
821     return m_faceTagsInUse.size();
822     default: {
823     stringstream msg;
824 caltinay 3791 msg << "getNumberOfTagsInUse: invalid function space type "
825     << fsType;
826 caltinay 5984 throw ValueError(msg.str());
827 caltinay 3722 }
828     }
829     }
830    
831     const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
832     {
833     switch(fsType) {
834     case Nodes:
835     return &m_nodeTagsInUse[0];
836     case Elements:
837     case ReducedElements:
838     return &m_elementTagsInUse[0];
839     case FaceElements:
840     case ReducedFaceElements:
841     return &m_faceTagsInUse[0];
842     default: {
843     stringstream msg;
844 caltinay 3791 msg << "borrowListOfTagsInUse: invalid function space type "
845     << fsType;
846 caltinay 5984 throw ValueError(msg.str());
847 caltinay 3722 }
848     }
849     }
850    
851 caltinay 5120 void RipleyDomain::Print_Mesh_Info(bool full) const
852 caltinay 3670 {
853 caltinay 3691 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
854     << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
855     cout << "Number of dimensions: " << m_numDim << endl;
856 jfenwick 4388 cout << "Number of elements per rank: " << getNumElements() << endl;
857 caltinay 3691 // write tags
858     if (m_tagMap.size() > 0) {
859     cout << "Tags:" << endl;
860     TagMap::const_iterator it;
861     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
862 caltinay 5120 cout << " " << setw(5) << it->second << " " << it->first << endl;
863 caltinay 3691 }
864 caltinay 3670 }
865     }
866    
867 caltinay 5148 int RipleyDomain::getSystemMatrixTypeId(const bp::object& options) const
868 caltinay 3670 {
869 caltinay 5148 const escript::SolverBuddy& sb = bp::extract<escript::SolverBuddy>(options);
870     int package = sb.getPackage();
871 caltinay 6395 escript::SolverOptions method = sb.getSolverMethod();
872 caltinay 6399 #ifdef ESYS_HAVE_TRILINOS
873 caltinay 6395 bool isDirect = escript::isDirectSolver(method);
874 caltinay 6399 #endif
875 caltinay 5148
876 caltinay 6166 // the configuration of ripley should have taken care that we have either
877     // paso or trilinos so here's how we prioritize
878 caltinay 6317 #if defined(ESYS_HAVE_PASO) && defined(ESYS_HAVE_TRILINOS)
879     // we have Paso & Trilinos so use Trilinos for parallel direct solvers
880 caltinay 6395 if (package == escript::SO_DEFAULT) {
881     if ((method == escript::SO_METHOD_DIRECT && getMPISize() > 1)
882     || isDirect
883     || sb.isComplex()) {
884     package = escript::SO_PACKAGE_TRILINOS;
885     }
886 caltinay 6317 }
887     #endif
888 caltinay 6166 #ifdef ESYS_HAVE_PASO
889     if (package == escript::SO_DEFAULT)
890 caltinay 5159 package = escript::SO_PACKAGE_PASO;
891     #endif
892 caltinay 6166 #ifdef ESYS_HAVE_TRILINOS
893     if (package == escript::SO_DEFAULT)
894     package = escript::SO_PACKAGE_TRILINOS;
895     #endif
896 aellery 6799 if (package == escript::SO_PACKAGE_TRILINOS) {
897 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
898 caltinay 6338 int type = (int)SMT_TRILINOS;
899     if (sb.isComplex())
900     type |= (int)SMT_COMPLEX;
901 caltinay 6395 // This is required because MueLu (AMG) and Amesos2 (direct) do not
902     // support block matrices at this point. Remove if they ever do...
903 caltinay 6376 if (sb.getPreconditioner() == escript::SO_PRECONDITIONER_AMG ||
904 caltinay 6377 sb.getPreconditioner() == escript::SO_PRECONDITIONER_ILUT ||
905 caltinay 6395 isDirect) {
906 caltinay 6376 type |= (int)SMT_UNROLL;
907     }
908 caltinay 6338 return type;
909 caltinay 6119 #else
910     throw RipleyException("Trilinos requested but not built with Trilinos.");
911     #endif
912 caltinay 5148 }
913 caltinay 6166 #ifdef ESYS_HAVE_PASO
914 caltinay 6338 if (sb.isComplex()) {
915     throw NotImplementedError("Paso does not support complex-valued matrices");
916     }
917 caltinay 5148 // in all other cases we use PASO
918 caltinay 5929 return (int)SMT_PASO | paso::SystemMatrix::getSystemMatrixTypeId(
919 caltinay 6395 method, sb.getPreconditioner(), sb.getPackage(),
920 caltinay 5148 sb.isSymmetric(), m_mpiInfo);
921 caltinay 6166 #else
922     throw RipleyException("Unable to find a working solver library!");
923     #endif
924 caltinay 3670 }
925    
926 caltinay 5120 int RipleyDomain::getTransportTypeId(int solver, int preconditioner,
927     int package, bool symmetry) const
928 caltinay 3670 {
929 caltinay 6166 #ifdef ESYS_USE_PASO
930     return paso::TransportProblem::getTypeId(solver, preconditioner,
931 caltinay 3691 package, symmetry, m_mpiInfo);
932 caltinay 6166 #else
933     throw RipleyException("Transport solvers require Paso but ripley was not compiled with Paso!");
934     #endif
935 caltinay 3670 }
936    
937 caltinay 5117 escript::ASM_ptr RipleyDomain::newSystemMatrix(int row_blocksize,
938     const escript::FunctionSpace& row_functionspace, int column_blocksize,
939     const escript::FunctionSpace& column_functionspace, int type) const
940 caltinay 3670 {
941 caltinay 3691 bool reduceRowOrder=false;
942     bool reduceColOrder=false;
943     // is the domain right?
944     const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
945 caltinay 3764 if (row_domain != *this)
946 caltinay 5984 throw ValueError("newSystemMatrix: domain of row function space does not match the domain of matrix generator");
947 caltinay 3691 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
948 caltinay 3764 if (col_domain != *this)
949 caltinay 5984 throw ValueError("newSystemMatrix: domain of column function space does not match the domain of matrix generator");
950 caltinay 3691 // is the function space type right?
951 caltinay 3750 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
952 caltinay 3691 reduceRowOrder=true;
953 caltinay 3750 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
954 caltinay 5984 throw ValueError("newSystemMatrix: illegal function space type for system matrix rows");
955 caltinay 3750 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
956 caltinay 3691 reduceColOrder=true;
957 caltinay 3750 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
958 caltinay 5984 throw ValueError("newSystemMatrix: illegal function space type for system matrix columns");
959 caltinay 5148 // are block sizes identical?
960     if (row_blocksize != column_blocksize)
961 caltinay 5984 throw ValueError("newSystemMatrix: row/column block sizes must be equal");
962 caltinay 5148 // are function spaces equal
963     if (reduceRowOrder != reduceColOrder)
964 caltinay 5984 throw ValueError("newSystemMatrix: row/column function spaces must be equal");
965 caltinay 3691
966 caltinay 5148 // generate matrix
967     //if (reduceRowOrder || reduceColOrder)
968 caltinay 5984 // throw NotImplementedError("newSystemMatrix: reduced order not supported");
969 caltinay 4775
970 caltinay 5209 if (type & (int)SMT_CUSP) {
971 aellery 6799 #ifndef ESYS_HAVE_CUDA
972     throw RipleyException("eScript does not support CUDA.");
973 caltinay 5159 #endif
974 caltinay 6119 } else if (type & (int)SMT_TRILINOS) {
975 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
976 caltinay 6119 const_TrilinosGraph_ptr graph(getTrilinosGraph());
977 caltinay 6338 bool isComplex = (type & (int)SMT_COMPLEX);
978 caltinay 6376 bool unroll = (type & (int)SMT_UNROLL);
979 caltinay 6119 escript::ASM_ptr sm(new TrilinosMatrixAdapter(m_mpiInfo, row_blocksize,
980 caltinay 6376 row_functionspace, graph, isComplex, unroll));
981 caltinay 6119 return sm;
982     #else
983     throw RipleyException("newSystemMatrix: ripley was not compiled with "
984     "Trilinos support so the Trilinos solver stack cannot be used.");
985     #endif
986 caltinay 5148 } else if (type & (int)SMT_PASO) {
987 caltinay 6166 #ifdef ESYS_HAVE_PASO
988 caltinay 5148 paso::SystemMatrixPattern_ptr pattern(getPasoMatrixPattern(
989     reduceRowOrder, reduceColOrder));
990     type -= (int)SMT_PASO;
991 caltinay 5929 escript::ASM_ptr sm(new paso::SystemMatrix(type, pattern,
992     row_blocksize, column_blocksize, false, row_functionspace,
993     column_functionspace));
994 caltinay 5148 return sm;
995 caltinay 6166 #else
996     throw RipleyException("newSystemMatrix: ripley was not compiled with "
997     "Paso support so the Paso solver stack cannot be used.");
998     #endif
999 caltinay 5148 } else {
1000     throw RipleyException("newSystemMatrix: unknown matrix type ID");
1001     }
1002 caltinay 3670 }
1003    
1004 caltinay 5114 void RipleyDomain::addToSystem(escript::AbstractSystemMatrix& mat,
1005     escript::Data& rhs, const DataMap& coefs,
1006     Assembler_ptr assembler) const
1007 sshaw 4622 {
1008     if (isNotEmpty("d_contact", coefs) || isNotEmpty("y_contact", coefs))
1009 caltinay 5984 throw ValueError(
1010 sshaw 4622 "addToSystem: Ripley does not support contact elements");
1011    
1012 caltinay 5114 assemblePDE(&mat, rhs, coefs, assembler);
1013     assemblePDEBoundary(&mat, rhs, coefs, assembler);
1014     assemblePDEDirac(&mat, rhs, coefs, assembler);
1015 sshaw 4622 }
1016    
1017     void RipleyDomain::addToSystemFromPython(escript::AbstractSystemMatrix& mat,
1018 caltinay 5120 escript::Data& rhs,
1019     const bp::list& data,
1020     Assembler_ptr assembler) const
1021 sshaw 4622 {
1022 caltinay 5114 DataMap mapping;
1023 sshaw 4622 tupleListToMap(mapping, data);
1024 sshaw 4942 addToSystem(mat, rhs, mapping, assembler);
1025 sshaw 4622 }
1026    
1027 caltinay 5120 Assembler_ptr RipleyDomain::createAssemblerFromPython(const string type,
1028     const bp::list& options) const
1029     {
1030 caltinay 5114 DataMap mapping;
1031 sshaw 4622 tupleListToMap(mapping, options);
1032 sshaw 4942 return createAssembler(type, mapping);
1033 caltinay 3760 }
1034    
1035 caltinay 5120 void RipleyDomain::addToRHSFromPython(escript::Data& rhs, const bp::list& data,
1036     Assembler_ptr assembler) const
1037 sshaw 4622 {
1038 caltinay 5114 DataMap mapping;
1039 sshaw 4622 tupleListToMap(mapping, data);
1040 sshaw 4942 addToRHS(rhs, mapping, assembler);
1041 sshaw 4622 }
1042    
1043 caltinay 5114 void RipleyDomain::addToRHS(escript::Data& rhs, const DataMap& coefs,
1044     Assembler_ptr assembler) const
1045 sshaw 4622 {
1046     if (isNotEmpty("y_contact", coefs))
1047 caltinay 5984 throw ValueError(
1048 sshaw 4622 "addPDEToRHS: Ripley does not support contact elements");
1049    
1050     if (rhs.isEmpty()) {
1051 sshaw 5573 if ((isNotEmpty("X", coefs) && isNotEmpty("du", coefs))
1052     || isNotEmpty("Y", coefs))
1053 caltinay 5984 throw ValueError(
1054 sshaw 4622 "addPDEToRHS: right hand side coefficients are provided "
1055     "but no right hand side vector given");
1056     else
1057     return;
1058     }
1059    
1060 caltinay 5114 assemblePDE(NULL, rhs, coefs, assembler);
1061     assemblePDEBoundary(NULL, rhs, coefs, assembler);
1062     assemblePDEDirac(NULL, rhs, coefs, assembler);
1063 sshaw 4622 }
1064    
1065 caltinay 5120 escript::ATP_ptr RipleyDomain::newTransportProblem(int blocksize,
1066     const escript::FunctionSpace& functionspace, int type) const
1067 caltinay 3836 {
1068     // is the domain right?
1069     const RipleyDomain& domain=dynamic_cast<const RipleyDomain&>(*(functionspace.getDomain()));
1070     if (domain != *this)
1071 caltinay 5984 throw ValueError("newTransportProblem: domain of function space does not match the domain of transport problem generator");
1072 caltinay 3836 // is the function space type right?
1073 caltinay 6169 if (functionspace.getTypeCode() != ReducedDegreesOfFreedom &&
1074     functionspace.getTypeCode() != DegreesOfFreedom)
1075 caltinay 5984 throw ValueError("newTransportProblem: illegal function space type for transport problem");
1076 caltinay 3836
1077 caltinay 6166 #ifdef ESYS_HAVE_PASO
1078 caltinay 6169 const bool reduced = (functionspace.getTypeCode() == ReducedDegreesOfFreedom);
1079 caltinay 3836 // generate matrix
1080 caltinay 6169 paso::SystemMatrixPattern_ptr pattern(getPasoMatrixPattern(reduced, reduced));
1081 caltinay 5929 escript::ATP_ptr tp(new paso::TransportProblem(pattern, blocksize,
1082     functionspace));
1083     return tp;
1084 caltinay 6166 #else
1085     throw RipleyException("newTransportProblem: transport problems require the "
1086     "Paso library which is not available.");
1087     #endif
1088 caltinay 3836 }
1089    
1090 sshaw 4942 void RipleyDomain::addPDEToTransportProblemFromPython(
1091 caltinay 5120 escript::AbstractTransportProblem& tp, escript::Data& source,
1092     const bp::list& data, Assembler_ptr assembler) const
1093 sshaw 4942 {
1094 caltinay 5114 DataMap mapping;
1095 sshaw 4942 tupleListToMap(mapping, data);
1096     addPDEToTransportProblem(tp, source, mapping, assembler);
1097     }
1098    
1099 caltinay 3836 void RipleyDomain::addPDEToTransportProblem(
1100 caltinay 5114 escript::AbstractTransportProblem& tp, escript::Data& source,
1101     const DataMap& coefs, Assembler_ptr assembler) const
1102 caltinay 3836 {
1103 sshaw 4942 if (isNotEmpty("d_contact", coefs) || isNotEmpty("y_contact", coefs))
1104 caltinay 5984 throw ValueError("addPDEToTransportProblem: Ripley does not support contact elements");
1105 caltinay 3836
1106 caltinay 6166 #ifdef ESYS_HAVE_PASO
1107     paso::TransportProblem* ptp = dynamic_cast<paso::TransportProblem*>(&tp);
1108 caltinay 5929 if (!ptp)
1109 caltinay 5984 throw ValueError("addPDEToTransportProblem: Ripley only accepts Paso transport problems");
1110 caltinay 3836
1111 caltinay 5929 escript::ASM_ptr mm(ptp->borrowMassMatrix());
1112     escript::ASM_ptr tm(ptp->borrowTransportMatrix());
1113 caltinay 5114
1114 caltinay 5929 assemblePDE(mm.get(), source, coefs, assembler);
1115     assemblePDE(tm.get(), source, coefs, assembler);
1116     assemblePDEBoundary(tm.get(), source, coefs, assembler);
1117     assemblePDEDirac(tm.get(), source, coefs, assembler);
1118 caltinay 6166 #else
1119     throw RipleyException("addPDEToTransportProblem: transport problems "
1120     "require the Paso library which is not available.");
1121     #endif
1122 caltinay 3836 }
1123    
1124 jfenwick 5303 void RipleyDomain::addPDEToTransportProblem(
1125     escript::AbstractTransportProblem& tp, escript::Data& source,
1126     const escript::Data& M,
1127     const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,
1128     const escript::Data& X,const escript::Data& Y,
1129     const escript::Data& d, const escript::Data& y,
1130     const escript::Data& d_contact,const escript::Data& y_contact,
1131     const escript::Data& d_dirac,const escript::Data& y_dirac) const
1132     {
1133     throw RipleyException("Programmer error: incorrect version of addPDEToTransportProblem called");
1134    
1135     }
1136    
1137 caltinay 3697 void RipleyDomain::setNewX(const escript::Data& arg)
1138     {
1139 caltinay 5984 throw NotImplementedError("setNewX(): operation not supported");
1140 caltinay 3697 }
1141    
1142 caltinay 3722 //protected
1143 caltinay 6580 template<typename Scalar>
1144 caltinay 6170 void RipleyDomain::dofToNodes(escript::Data& out, const escript::Data& in) const
1145     {
1146     // expand data object if necessary
1147     const_cast<escript::Data*>(&in)->expand();
1148     const dim_t numComp = in.getDataPointSize();
1149     const dim_t numNodes = getNumNodes();
1150 caltinay 6581 const Scalar zero = static_cast<Scalar>(0);
1151 caltinay 6170 out.requireWrite();
1152 caltinay 6580 #ifdef ESYS_HAVE_TRILINOS
1153 caltinay 6170 using namespace esys_trilinos;
1154    
1155 caltinay 6213 const_TrilinosGraph_ptr graph(getTrilinosGraph());
1156     Teuchos::RCP<const MapType> colMap;
1157     Teuchos::RCP<const MapType> rowMap;
1158     MapType colPointMap;
1159     MapType rowPointMap;
1160 caltinay 6170
1161 caltinay 6213 if (numComp > 1) {
1162 caltinay 6580 colPointMap = BlockVectorType<Scalar>::makePointMap(
1163     *graph->getColMap(), numComp);
1164     rowPointMap = BlockVectorType<Scalar>::makePointMap(
1165     *graph->getRowMap(), numComp);
1166 caltinay 6213 colMap = Teuchos::rcpFromRef(colPointMap);
1167     rowMap = Teuchos::rcpFromRef(rowPointMap);
1168     } else {
1169     colMap = graph->getColMap();
1170     rowMap = graph->getRowMap();
1171     }
1172    
1173     const ImportType importer(rowMap, colMap);
1174 caltinay 6580 const Teuchos::ArrayView<const Scalar> localIn(
1175     in.getSampleDataRO(0, zero), in.getNumDataPoints()*numComp);
1176     Teuchos::RCP<VectorType<Scalar> > lclData = rcp(new VectorType<Scalar>(
1177     rowMap, localIn, localIn.size(), 1));
1178     Teuchos::RCP<VectorType<Scalar> > gblData = rcp(new VectorType<Scalar>(
1179     colMap, 1));
1180 caltinay 6170 gblData->doImport(*lclData, importer, Tpetra::INSERT);
1181 caltinay 6580 Teuchos::ArrayRCP<const Scalar> gblArray(gblData->getData(0));
1182 caltinay 6213
1183 caltinay 6170 #pragma omp parallel for
1184     for (index_t i = 0; i < numNodes; i++) {
1185 caltinay 6580 const Scalar* src = &gblArray[getDofOfNode(i) * numComp];
1186     copy(src, src+numComp, out.getSampleDataRW(i, zero));
1187     }
1188     #elif defined(ESYS_HAVE_PASO)
1189 caltinay 6581 paso::Coupler_ptr<Scalar> coupler(new paso::Coupler<Scalar>(m_connector,
1190     numComp, m_mpiInfo));
1191     coupler->startCollect(in.getSampleDataRO(0, zero));
1192 caltinay 6580 const dim_t numDOF = getNumDOF();
1193 caltinay 6581 const Scalar* buffer = coupler->finishCollect();
1194 caltinay 6580
1195     #pragma omp parallel for
1196     for (index_t i = 0; i < numNodes; i++) {
1197     const index_t dof = getDofOfNode(i);
1198 caltinay 6581 const Scalar* src = (dof < numDOF ? in.getSampleDataRO(dof, zero)
1199 caltinay 6580 : &buffer[(dof - numDOF) * numComp]);
1200 caltinay 6581 copy(src, src+numComp, out.getSampleDataRW(i, zero));
1201 caltinay 6170 }
1202 caltinay 6580 #endif // ESYS_HAVE_PASO
1203 caltinay 6170 }
1204    
1205     //protected
1206 caltinay 6435 template<typename Scalar>
1207 caltinay 4626 void RipleyDomain::copyData(escript::Data& out, const escript::Data& in) const
1208 caltinay 4211 {
1209     const dim_t numComp = in.getDataPointSize();
1210 caltinay 5039 const dim_t numSamples = in.getNumSamples();
1211 caltinay 6435 const Scalar zero = static_cast<Scalar>(0);
1212 caltinay 4211 out.requireWrite();
1213     #pragma omp parallel for
1214 caltinay 6435 for (index_t i = 0; i < numSamples; i++) {
1215     const Scalar* src = in.getSampleDataRO(i, zero);
1216     copy(src, src+numComp, out.getSampleDataRW(i, zero));
1217 caltinay 4211 }
1218     }
1219    
1220     //protected
1221 caltinay 6435 template<typename Scalar>
1222 caltinay 4626 void RipleyDomain::averageData(escript::Data& out, const escript::Data& in) const
1223 caltinay 3722 {
1224     const dim_t numComp = in.getDataPointSize();
1225 caltinay 3744 const dim_t dpp = in.getNumDataPointsPerSample();
1226 caltinay 5039 const dim_t numSamples = in.getNumSamples();
1227 caltinay 6435 const Scalar zero = static_cast<Scalar>(0);
1228 caltinay 3740 out.requireWrite();
1229 caltinay 3722 #pragma omp parallel for
1230 caltinay 6435 for (index_t i = 0; i < numSamples; i++) {
1231     const Scalar* src = in.getSampleDataRO(i, zero);
1232     Scalar* dest = out.getSampleDataRW(i, zero);
1233     for (index_t c = 0; c < numComp; c++) {
1234     Scalar res = zero;
1235     for (index_t q = 0; q < dpp; q++)
1236     res += src[c+q*numComp];
1237     *dest++ = res / static_cast<real_t>(dpp);
1238 caltinay 3744 }
1239     }
1240     }
1241    
1242     //protected
1243 caltinay 6435 template<typename Scalar>
1244 caltinay 4626 void RipleyDomain::multiplyData(escript::Data& out, const escript::Data& in) const
1245 caltinay 3744 {
1246     const dim_t numComp = in.getDataPointSize();
1247 caltinay 4211 const dim_t dpp = out.getNumDataPointsPerSample();
1248 caltinay 5039 const dim_t numSamples = in.getNumSamples();
1249 caltinay 6435 const Scalar zero = static_cast<Scalar>(0);
1250 caltinay 3744 out.requireWrite();
1251     #pragma omp parallel for
1252 caltinay 6435 for (index_t i = 0; i < numSamples; i++) {
1253     const Scalar* src = in.getSampleDataRO(i, zero);
1254     Scalar* dest = out.getSampleDataRW(i, zero);
1255     for (index_t c = 0; c < numComp; c++) {
1256     for (index_t q = 0; q < dpp; q++)
1257 caltinay 4211 dest[c+q*numComp] = src[c];
1258     }
1259 caltinay 3722 }
1260     }
1261    
1262     //protected
1263     void RipleyDomain::updateTagsInUse(int fsType) const
1264     {
1265 caltinay 5187 vector<int>* tagsInUse=NULL;
1266     const vector<int>* tags=NULL;
1267 caltinay 3722 switch(fsType) {
1268     case Nodes:
1269     tags=&m_nodeTags;
1270     tagsInUse=&m_nodeTagsInUse;
1271     break;
1272     case Elements:
1273     case ReducedElements:
1274     tags=&m_elementTags;
1275     tagsInUse=&m_elementTagsInUse;
1276     break;
1277     case FaceElements:
1278     case ReducedFaceElements:
1279     tags=&m_faceTags;
1280     tagsInUse=&m_faceTagsInUse;
1281     break;
1282 sshaw 4660 case Points:
1283 caltinay 5984 throw NotImplementedError("updateTagsInUse for Ripley dirac points"
1284     " not supported");
1285 caltinay 3722 default:
1286     return;
1287     }
1288    
1289     // gather global unique tag values from tags into tagsInUse
1290     tagsInUse->clear();
1291 caltinay 5187 int lastFoundValue = numeric_limits<int>::min();
1292     int minFoundValue, local_minFoundValue;
1293     const int numTags = tags->size();
1294 caltinay 3722
1295     while (true) {
1296 caltinay 5120 // find smallest value bigger than lastFoundValue
1297 caltinay 5187 minFoundValue = numeric_limits<int>::max();
1298 caltinay 3722 #pragma omp parallel private(local_minFoundValue)
1299     {
1300     local_minFoundValue = minFoundValue;
1301 caltinay 5187 #pragma omp for schedule(static) nowait
1302     for (int i = 0; i < numTags; i++) {
1303     const int v = (*tags)[i];
1304 caltinay 3722 if ((v > lastFoundValue) && (v < local_minFoundValue))
1305     local_minFoundValue = v;
1306     }
1307     #pragma omp critical
1308     {
1309     if (local_minFoundValue < minFoundValue)
1310     minFoundValue = local_minFoundValue;
1311     }
1312     }
1313     #ifdef ESYS_MPI
1314     local_minFoundValue = minFoundValue;
1315     MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
1316     #endif
1317 caltinay 3740
1318 caltinay 3722 // if we found a new value add it to the tagsInUse vector
1319 caltinay 5187 if (minFoundValue < numeric_limits<int>::max()) {
1320 caltinay 3722 tagsInUse->push_back(minFoundValue);
1321     lastFoundValue = minFoundValue;
1322     } else
1323     break;
1324     }
1325     }
1326    
1327 caltinay 6166 #ifdef ESYS_HAVE_PASO
1328 caltinay 6170 void RipleyDomain::createPasoConnector(const RankVector& neighbour,
1329     const IndexVector& offsetInSharedSend,
1330     const IndexVector& offsetInSharedRecv,
1331     const IndexVector& sendShared,
1332     const IndexVector& recvShared)
1333     {
1334 caltinay 6172 const index_t* sendPtr = neighbour.empty() ? NULL : &sendShared[0];
1335     const index_t* recvPtr = neighbour.empty() ? NULL : &recvShared[0];
1336 caltinay 6170 paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
1337 caltinay 6173 getNumDOF(), neighbour, sendPtr, offsetInSharedSend));
1338 caltinay 6170 paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
1339 caltinay 6173 getNumDOF(), neighbour, recvPtr, offsetInSharedRecv));
1340 caltinay 6170 m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
1341     }
1342    
1343 caltinay 3756 //protected
1344 caltinay 5117 paso::Pattern_ptr RipleyDomain::createPasoPattern(
1345 caltinay 5120 const vector<IndexVector>& indices, dim_t N) const
1346 caltinay 3756 {
1347     // paso will manage the memory
1348 caltinay 5117 const dim_t M = indices.size();
1349     index_t* ptr = new index_t[M+1];
1350     ptr[0] = 0;
1351     for (index_t i=0; i < M; i++) {
1352     ptr[i+1] = ptr[i]+indices[i].size();
1353 caltinay 5114 }
1354    
1355 caltinay 5117 index_t* index = new index_t[ptr[M]];
1356 caltinay 3756
1357 caltinay 5117 #pragma omp parallel for
1358     for (index_t i=0; i < M; i++) {
1359     copy(indices[i].begin(), indices[i].end(), &index[ptr[i]]);
1360 caltinay 3756 }
1361    
1362 caltinay 5117 return paso::Pattern_ptr(new paso::Pattern(MATRIX_FORMAT_DEFAULT, M, N, ptr, index));
1363 caltinay 3756 }
1364 caltinay 6166 #endif // ESYS_HAVE_PASO
1365 caltinay 3756
1366 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
1367 caltinay 3756 //protected
1368 caltinay 6119 esys_trilinos::const_TrilinosGraph_ptr RipleyDomain::createTrilinosGraph(
1369     const IndexVector& myRows,
1370     const IndexVector& myColumns) const
1371     {
1372     using namespace esys_trilinos;
1373    
1374     const dim_t numMatrixRows = getNumDOF();
1375    
1376     TrilinosMap_ptr rowMap(new MapType(getNumDataPointsGlobal(), myRows,
1377     0, TeuchosCommFromEsysComm(m_mpiInfo->comm)));
1378    
1379     IndexVector columns(getNumNodes());
1380     // order is important - our columns (=myRows) come first, followed by
1381     // shared ones (=node Id for non-DOF)
1382     #pragma omp parallel for
1383     for (size_t i=0; i<columns.size(); i++) {
1384     columns[getDofOfNode(i)] = myColumns[i];
1385     }
1386     TrilinosMap_ptr colMap(new MapType(getNumDataPointsGlobal(), columns,
1387     0, TeuchosCommFromEsysComm(m_mpiInfo->comm)));
1388    
1389     // now build CSR arrays (rowPtr and colInd)
1390     const vector<IndexVector>& conns(getConnections(true));
1391     Teuchos::ArrayRCP<size_t> rowPtr(numMatrixRows+1);
1392     for (size_t i=0; i < numMatrixRows; i++) {
1393     rowPtr[i+1] = rowPtr[i] + conns[i].size();
1394     }
1395    
1396     Teuchos::ArrayRCP<LO> colInd(rowPtr[numMatrixRows]);
1397    
1398     #pragma omp parallel for
1399     for (index_t i=0; i < numMatrixRows; i++) {
1400     copy(conns[i].begin(), conns[i].end(), &colInd[rowPtr[i]]);
1401     }
1402    
1403     TrilinosGraph_ptr graph(new GraphType(rowMap, colMap, rowPtr, colInd));
1404     Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::parameterList();
1405     params->set("Optimize Storage", true);
1406     graph->fillComplete(rowMap, rowMap, params);
1407     return graph;
1408     }
1409 caltinay 6166 #endif // ESYS_HAVE_TRILINOS
1410 caltinay 6119
1411     //protected
1412 caltinay 6338 template<>
1413     void RipleyDomain::addToSystemMatrix<real_t>(escript::AbstractSystemMatrix* mat,
1414     const IndexVector& nodes, dim_t numEq,
1415     const DoubleVector& array) const
1416 caltinay 5114 {
1417 caltinay 6166 #ifdef ESYS_HAVE_PASO
1418 caltinay 6119 paso::SystemMatrix* psm = dynamic_cast<paso::SystemMatrix*>(mat);
1419     if (psm) {
1420     addToPasoMatrix(psm, nodes, numEq, array);
1421     return;
1422     }
1423 caltinay 6166 #endif
1424 caltinay 6144 #ifdef ESYS_HAVE_CUDA
1425 caltinay 6119 SystemMatrix* rsm = dynamic_cast<SystemMatrix*>(mat);
1426     if (rsm) {
1427     rsm->add(nodes, array);
1428     return;
1429     }
1430 caltinay 5162 #endif
1431 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
1432 caltinay 6119 TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(mat);
1433     if (tm) {
1434     tm->add(nodes, array);
1435     return;
1436 caltinay 5114 }
1437 caltinay 6119 #endif
1438     throw RipleyException("addToSystemMatrix: unknown system matrix type");
1439 caltinay 5114 }
1440    
1441 caltinay 6338 //protected
1442     template<>
1443     void RipleyDomain::addToSystemMatrix<cplx_t>(escript::AbstractSystemMatrix* mat,
1444     const IndexVector& nodes, dim_t numEq,
1445     const vector<cplx_t>& array) const
1446     {
1447     #ifdef ESYS_HAVE_TRILINOS
1448     TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(mat);
1449     if (tm) {
1450     tm->add(nodes, array);
1451     return;
1452     }
1453     #endif
1454     throw RipleyException("addToSystemMatrix: only Trilinos matrices support "
1455     "complex-valued assembly!");
1456     }
1457    
1458 caltinay 6166 #ifdef ESYS_HAVE_PASO
1459 caltinay 5114 //private
1460 caltinay 5929 void RipleyDomain::addToPasoMatrix(paso::SystemMatrix* mat,
1461     const IndexVector& nodes, dim_t numEq,
1462     const vector<double>& array) const
1463 caltinay 3756 {
1464     const dim_t numMyCols = mat->pattern->mainPattern->numInput;
1465     const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
1466 caltinay 5114 const dim_t numSubblocks_Eq = numEq / mat->row_block_size;
1467     const dim_t numSubblocks_Sol = numEq / mat->col_block_size;
1468 caltinay 3756
1469     const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
1470     const index_t* mainBlock_index = mat->mainBlock->pattern->index;
1471     double* mainBlock_val = mat->mainBlock->val;
1472     const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
1473     const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
1474     double* col_coupleBlock_val = mat->col_coupleBlock->val;
1475     const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
1476     const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
1477     double* row_coupleBlock_val = mat->row_coupleBlock->val;
1478 caltinay 3841 index_t offset=(mat->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
1479 caltinay 3756
1480 caltinay 3850 #define UPDATE_BLOCK(VAL) do {\
1481     for (dim_t ic=0; ic<mat->col_block_size; ++ic) {\
1482     const dim_t i_Sol=ic+mat->col_block_size*l_col;\
1483     for (dim_t ir=0; ir<mat->row_block_size; ++ir) {\
1484     const dim_t i_Eq=ir+mat->row_block_size*l_row;\
1485     VAL[k*mat->block_size+ir+mat->row_block_size*ic]\
1486 caltinay 5114 += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, numEq, numEq, nodes.size())];\
1487 caltinay 3850 }\
1488     }\
1489     } while(0)
1490    
1491     if (mat->type & MATRIX_FORMAT_CSC) {
1492 caltinay 5114 for (dim_t k_Sol = 0; k_Sol < nodes.size(); ++k_Sol) {
1493 caltinay 3850 // down columns of array
1494     for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
1495 caltinay 5114 const dim_t i_col = nodes[k_Sol]*numSubblocks_Sol+l_col;
1496 caltinay 3850 if (i_col < numMyCols) {
1497 caltinay 5114 for (dim_t k_Eq = 0; k_Eq < nodes.size(); ++k_Eq) {
1498 caltinay 3850 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
1499 caltinay 5114 const dim_t i_row = nodes[k_Eq]*numSubblocks_Eq+l_row+offset;
1500 caltinay 3850 if (i_row < numMyRows+offset) {
1501     for (dim_t k = mainBlock_ptr[i_col]-offset; k < mainBlock_ptr[i_col+1]-offset; ++k) {
1502     if (mainBlock_index[k] == i_row) {
1503     UPDATE_BLOCK(mainBlock_val);
1504     break;
1505 caltinay 3760 }
1506     }
1507 caltinay 3850 } else {
1508     for (dim_t k = col_coupleBlock_ptr[i_col]-offset; k < col_coupleBlock_ptr[i_col+1]-offset; ++k) {
1509     if (row_coupleBlock_index[k] == i_row - numMyRows) {
1510     UPDATE_BLOCK(row_coupleBlock_val);
1511     break;
1512     }
1513     }
1514 caltinay 3760 }
1515 caltinay 3850 }
1516     }
1517     } else {
1518 caltinay 5114 for (dim_t k_Eq = 0; k_Eq < nodes.size(); ++k_Eq) {
1519 caltinay 3850 // across rows of array
1520     for (dim_t l_row=0; l_row<numSubblocks_Eq; ++l_row) {
1521 caltinay 5114 const dim_t i_row = nodes[k_Eq]*numSubblocks_Eq+l_row+offset;
1522 caltinay 3850 if (i_row < numMyRows+offset) {
1523     for (dim_t k = col_coupleBlock_ptr[i_col-numMyCols]-offset;
1524     k < col_coupleBlock_ptr[i_col-numMyCols+1]-offset; ++k)
1525     {
1526     if (col_coupleBlock_index[k] == i_row) {
1527     UPDATE_BLOCK(col_coupleBlock_val);
1528     break;
1529 caltinay 3760 }
1530     }
1531     }
1532 caltinay 3756 }
1533     }
1534     }
1535 caltinay 3850 }
1536     }
1537     } else {
1538 caltinay 5114 for (dim_t k_Eq = 0; k_Eq < nodes.size(); ++k_Eq) {
1539 caltinay 3850 // down columns of array
1540     for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
1541 caltinay 5114 const dim_t i_row = nodes[k_Eq]*numSubblocks_Eq+l_row;
1542 caltinay 3850 // only look at the matrix rows stored on this processor
1543     if (i_row < numMyRows) {
1544 caltinay 5114 for (dim_t k_Sol = 0; k_Sol < nodes.size(); ++k_Sol) {
1545 caltinay 3850 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
1546 caltinay 5114 const dim_t i_col = nodes[k_Sol]*numSubblocks_Sol+l_col+offset;
1547 caltinay 3850 if (i_col < numMyCols+offset) {
1548     for (dim_t k = mainBlock_ptr[i_row]-offset; k < mainBlock_ptr[i_row+1]-offset; ++k) {
1549     if (mainBlock_index[k] == i_col) {
1550     UPDATE_BLOCK(mainBlock_val);
1551     break;
1552 caltinay 3760 }
1553     }
1554 caltinay 3850 } else {
1555     for (dim_t k = col_coupleBlock_ptr[i_row]-offset; k < col_coupleBlock_ptr[i_row+1]-offset; ++k) {
1556     if (col_coupleBlock_index[k] == i_col-numMyCols) {
1557     UPDATE_BLOCK(col_coupleBlock_val);
1558     break;
1559     }
1560     }
1561 caltinay 3760 }
1562 caltinay 3756 }
1563     }
1564 caltinay 3850 } else {
1565 caltinay 5114 for (dim_t k_Sol = 0; k_Sol < nodes.size(); ++k_Sol) {
1566 caltinay 3850 // across rows of array
1567     for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1568 caltinay 5114 const dim_t i_col = nodes[k_Sol]*numSubblocks_Sol+l_col+offset;
1569 caltinay 3850 if (i_col < numMyCols+offset) {
1570     for (dim_t k = row_coupleBlock_ptr[i_row-numMyRows]-offset;
1571     k < row_coupleBlock_ptr[i_row-numMyRows+1]-offset; ++k)
1572     {
1573     if (row_coupleBlock_index[k] == i_col) {
1574     UPDATE_BLOCK(row_coupleBlock_val);
1575     break;
1576     }
1577     }
1578     }
1579     }
1580     }
1581 caltinay 3756 }
1582     }
1583     }
1584     }
1585 caltinay 3850 #undef UPDATE_BLOCK
1586 caltinay 3756 }
1587 caltinay 6166 #endif // ESYS_HAVE_PASO
1588 caltinay 3756
1589 caltinay 3836 //private
1590 caltinay 5114 void RipleyDomain::assemblePDE(escript::AbstractSystemMatrix* mat,
1591     escript::Data& rhs, const DataMap& coefs,
1592     Assembler_ptr assembler) const
1593 caltinay 3836 {
1594 sshaw 5573 if (rhs.isEmpty() && (isNotEmpty("X", coefs) || isNotEmpty("du", coefs))
1595     && isNotEmpty("Y", coefs))
1596 caltinay 5984 throw ValueError("assemblePDE: right hand side coefficients are "
1597     "provided but no right hand side vector given");
1598 caltinay 3836
1599     vector<int> fsTypes;
1600 sshaw 4760 assembler->collateFunctionSpaceTypes(fsTypes, coefs);
1601 caltinay 5120
1602 sshaw 4622 if (fsTypes.empty()) {
1603 caltinay 3836 return;
1604 sshaw 4622 }
1605 sshaw 5573
1606 caltinay 3836 int fs=fsTypes[0];
1607     if (fs != Elements && fs != ReducedElements)
1608 caltinay 5984 throw ValueError("assemblePDE: illegal function space type for coefficients");
1609 caltinay 3836
1610     for (vector<int>::const_iterator it=fsTypes.begin()+1; it!=fsTypes.end(); it++) {
1611     if (*it != fs) {
1612 caltinay 5984 throw ValueError("assemblePDE: coefficient function spaces don't match");
1613 caltinay 3836 }
1614     }
1615    
1616 caltinay 4023 int numEq, numComp;
1617     if (!mat) {
1618     if (rhs.isEmpty()) {
1619     numEq=numComp=1;
1620     } else {
1621     numEq=numComp=rhs.getDataPointSize();
1622     }
1623     } else {
1624 caltinay 5114 if (!rhs.isEmpty() && rhs.getDataPointSize() != mat->getRowBlockSize())
1625 caltinay 5984 throw ValueError("assemblePDE: matrix row block size and number of components of right hand side don't match");
1626 caltinay 5114 numEq = mat->getRowBlockSize();
1627     numComp = mat->getColumnBlockSize();
1628 caltinay 4023 }
1629 caltinay 3836
1630     if (numEq != numComp)
1631 caltinay 5984 throw ValueError("assemblePDE: number of equations and number of solutions don't match");
1632 caltinay 3836
1633 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
1634 caltinay 6119 TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(mat);
1635     if (tm) {
1636     tm->resumeFill();
1637     }
1638     #endif
1639 caltinay 3836
1640     if (numEq==1) {
1641 sshaw 4622 if (fs==ReducedElements) {
1642     assembler->assemblePDESingleReduced(mat, rhs, coefs);
1643     } else {
1644     assembler->assemblePDESingle(mat, rhs, coefs);
1645     }
1646 caltinay 3836 } else {
1647 sshaw 4622 if (fs==ReducedElements) {
1648     assembler->assemblePDESystemReduced(mat, rhs, coefs);
1649     } else {
1650     assembler->assemblePDESystem(mat, rhs, coefs);
1651     }
1652 caltinay 3836 }
1653 caltinay 6119
1654 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
1655 caltinay 6119 if (tm) {
1656     tm->fillComplete(true);
1657     }
1658     #endif
1659 caltinay 3836 }
1660    
1661     //private
1662 caltinay 5114 void RipleyDomain::assemblePDEBoundary(escript::AbstractSystemMatrix* mat,
1663     escript::Data& rhs, const DataMap& coefs, Assembler_ptr assembler) const
1664 caltinay 3836 {
1665 sshaw 4622 if (rhs.isEmpty() && isNotEmpty("y", coefs))
1666 caltinay 5984 throw ValueError("assemblePDEBoundary: y provided but no right hand side vector given");
1667 caltinay 3836
1668     int fs=-1;
1669 sshaw 4622 if (isNotEmpty("d", coefs))
1670 caltinay 5114 fs=coefs.find("d")->second.getFunctionSpace().getTypeCode();
1671 sshaw 4622 if (isNotEmpty("y", coefs)) {
1672 caltinay 5114 DataMap::const_iterator iy = coefs.find("y");
1673 caltinay 3836 if (fs == -1)
1674 sshaw 4622 fs = iy->second.getFunctionSpace().getTypeCode();
1675     else if (fs != iy->second.getFunctionSpace().getTypeCode())
1676 caltinay 5984 throw ValueError("assemblePDEBoundary: coefficient function spaces don't match");
1677 caltinay 3836 }
1678 sshaw 4622 if (fs==-1) {
1679     return;
1680     }
1681 caltinay 5120
1682 caltinay 3836 if (fs != FaceElements && fs != ReducedFaceElements)
1683 caltinay 5984 throw ValueError("assemblePDEBoundary: illegal function space type for coefficients");
1684 caltinay 3836
1685 caltinay 4648 int numEq, numComp;
1686     if (!mat) {
1687     if (rhs.isEmpty()) {
1688     numEq=numComp=1;
1689     } else {
1690     numEq=numComp=rhs.getDataPointSize();
1691     }
1692     } else {
1693 caltinay 5114 if (!rhs.isEmpty() && rhs.getDataPointSize() != mat->getRowBlockSize())
1694 caltinay 5984 throw ValueError("assemblePDEBoundary: matrix row block size and number of components of right hand side don't match");
1695 caltinay 5114 numEq = mat->getRowBlockSize();
1696     numComp = mat->getColumnBlockSize();
1697 caltinay 4648 }
1698 caltinay 3836
1699     if (numEq != numComp)
1700 caltinay 5984 throw ValueError("assemblePDEBoundary: number of equations and number of solutions don't match");
1701 caltinay 3836
1702 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
1703 caltinay 6119 TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(mat);
1704     if (tm) {
1705     tm->resumeFill();
1706     }
1707     #endif
1708 caltinay 3836
1709     if (numEq==1) {
1710     if (fs==ReducedFaceElements)
1711 sshaw 4622 assembler->assemblePDEBoundarySingleReduced(mat, rhs, coefs);
1712 caltinay 3836 else
1713 sshaw 4622 assembler->assemblePDEBoundarySingle(mat, rhs, coefs);
1714 caltinay 3836 } else {
1715     if (fs==ReducedFaceElements)
1716 sshaw 4622 assembler->assemblePDEBoundarySystemReduced(mat, rhs, coefs);
1717 caltinay 3836 else
1718 sshaw 4622 assembler->assemblePDEBoundarySystem(mat, rhs, coefs);
1719 caltinay 3836 }
1720 caltinay 6119
1721 caltinay 6144 #ifdef ESYS_HAVE_TRILINOS
1722 caltinay 6119 if (tm) {
1723     tm->fillComplete(true);
1724     }
1725     #endif
1726 caltinay 3836 }
1727    
1728 caltinay 5114 void RipleyDomain::assemblePDEDirac(escript::AbstractSystemMatrix* mat,
1729     escript::Data& rhs, const DataMap& coefs,
1730     Assembler_ptr assembler) const
1731 sshaw 4622 {
1732 caltinay 5114 bool yNotEmpty = isNotEmpty("y_dirac", coefs);
1733     bool dNotEmpty = isNotEmpty("d_dirac", coefs);
1734 sshaw 4622 if (!(yNotEmpty || dNotEmpty)) {
1735     return;
1736     }
1737 caltinay 5114 escript::Data d = unpackData("d_dirac", coefs);
1738     escript::Data y = unpackData("y_dirac", coefs);
1739 sshaw 4622 int nEq, nComp;
1740     if (!mat) {
1741     if (rhs.isEmpty()) {
1742     nEq=nComp=1;
1743     } else {
1744     nEq=nComp=rhs.getDataPointSize();
1745     }
1746     } else {
1747 caltinay 5114 if (!rhs.isEmpty() && rhs.getDataPointSize() != mat->getRowBlockSize())
1748 caltinay 5984 throw ValueError("assemblePDEDirac: matrix row block size "
1749 sshaw 4622 "and number of components of right hand side don't match");
1750 caltinay 5114 nEq = mat->getRowBlockSize();
1751     nComp = mat->getColumnBlockSize();
1752 sshaw 4622 }
1753 caltinay 5114
1754 sshaw 5582 rhs.requireWrite();
1755 sshaw 4629 for (int i = 0; i < m_diracPoints.size(); i++) { //only for this rank
1756 caltinay 5254 const index_t dof = getDofOfNode(m_diracPoints[i].node);
1757 sshaw 4622 if (yNotEmpty) {
1758 caltinay 4626 const double *EM_F = y.getSampleDataRO(i);
1759     double *F_p = rhs.getSampleDataRW(0);
1760 caltinay 5254 if (dof < getNumDOF()) {
1761 caltinay 4807 for (index_t eq = 0; eq < nEq; eq++) {
1762 caltinay 5254 F_p[INDEX2(eq, dof, nEq)] += EM_F[eq];
1763 caltinay 4807 }
1764 sshaw 4622 }
1765     }
1766     if (dNotEmpty) {
1767 caltinay 5254 const IndexVector rowIndex(1, dof);
1768 caltinay 4626 const double *EM_S = d.getSampleDataRO(i);
1769 caltinay 5120 vector<double> contents(EM_S, EM_S+nEq*nEq*nComp);
1770 caltinay 5114 addToSystemMatrix(mat, rowIndex, nEq, contents);
1771 sshaw 4622 }
1772     }
1773     }
1774    
1775 sshaw 5191 bool RipleyDomain::probeInterpolationAcross(int fsType_source,
1776 caltinay 5120 const escript::AbstractDomain&, int fsType_target) const
1777 caltinay 3670 {
1778 caltinay 4284 return false;
1779 caltinay 3670 }
1780    
1781 sshaw 5191 void RipleyDomain::interpolateAcross(escript::Data& target,
1782 caltinay 5120 const escript::Data& source) const
1783 caltinay 3670 {
1784 caltinay 5984 throw NotImplementedError("interpolateAcross() not supported");
1785 caltinay 3670 }
1786    
1787 jfenwick 4521 // Expecting ("gaussian", radius, sigma)
1788 caltinay 5120 bool RipleyDomain::supportsFilter(const bp::tuple& t) const
1789 jfenwick 4521 {
1790 caltinay 5120 if (len(t) == 0) { // so we can handle unfiltered randoms
1791 jfenwick 4687 return true;
1792     }
1793 caltinay 5120 if (len(t) != 3) {
1794 jfenwick 4521 return false;
1795     }
1796 caltinay 5120 bp::extract<string> ex(t[0]);
1797     if (!ex.check() || (ex() != "gaussian")) {
1798 jfenwick 4521 return false;
1799     }
1800 caltinay 5120 if (! bp::extract<unsigned int>(t[1]).check()) {
1801 jfenwick 4521 return false;
1802     }
1803 caltinay 5120 return bp::extract<double>(t[2]).check();
1804 jfenwick 4521 }
1805 caltinay 3670
1806 caltinay 5120 void RipleyDomain::addPoints(const vector<double>& coords,
1807     const vector<int>& tags)
1808 sshaw 4622 {
1809 caltinay 5097 for (int i = 0; i < tags.size(); i++) {
1810 caltinay 5120 dim_t node = findNode(&coords[i * m_numDim]);
1811 sshaw 4622 if (node >= 0) {
1812 sshaw 4660 m_diracPointNodeIDs.push_back(borrowSampleReferenceIDs(Nodes)[node]);
1813 caltinay 4626 DiracPoint dp;
1814 sshaw 4660 dp.node = node; //local
1815 caltinay 5097 dp.tag = tags[i];
1816 sshaw 4629 m_diracPoints.push_back(dp);
1817 sshaw 4622 }
1818     }
1819     }
1820 jfenwick 4521
1821 caltinay 3670 } // end of namespace ripley
1822    

  ViewVC Help
Powered by ViewVC 1.1.26