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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3766 - (hide annotations)
Thu Jan 12 08:17:49 2012 UTC (7 years, 8 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 45120 byte(s)
Checkpoint. Getting weipa to like ripley again.

1 caltinay 3670
2     /*******************************************************
3     *
4 caltinay 3764 * Copyright (c) 2003-2012 by University of Queensland
5 caltinay 3670 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14     #include <ripley/RipleyDomain.h>
15     #include <escript/DataFactory.h>
16 caltinay 3691 #include <escript/FunctionSpaceFactory.h>
17     #include <pasowrap/SystemMatrixAdapter.h>
18     #include <pasowrap/TransportProblemAdapter.h>
19 caltinay 3670
20     #include <iomanip>
21    
22     using namespace std;
23 caltinay 3691 using paso::SystemMatrixAdapter;
24     using paso::TransportProblemAdapter;
25 caltinay 3670
26     namespace ripley {
27    
28 caltinay 3691 escript::Domain_ptr RipleyDomain::loadMesh(const string& filename)
29 caltinay 3670 {
30 caltinay 3691 throw RipleyException("loadMesh() not implemented");
31 caltinay 3670 }
32    
33 caltinay 3691 escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34 caltinay 3670 {
35 caltinay 3691 throw RipleyException("readMesh() not implemented");
36 caltinay 3670 }
37    
38 caltinay 3691 RipleyDomain::RipleyDomain(dim_t dim) :
39     m_numDim(dim),
40     m_status(0)
41 caltinay 3670 {
42 caltinay 3691 m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
43 caltinay 3670 }
44    
45 caltinay 3691 RipleyDomain::~RipleyDomain()
46 caltinay 3670 {
47 caltinay 3691 Esys_MPIInfo_free(m_mpiInfo);
48 caltinay 3670 }
49    
50 caltinay 3744 bool RipleyDomain::operator==(const AbstractDomain& other) const
51     {
52     const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
53     if (o) {
54     return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
55     && m_elementTags==o->m_elementTags
56     && m_faceTags==o->m_faceTags);
57     }
58     return false;
59     }
60    
61 caltinay 3691 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
62 caltinay 3670 {
63 caltinay 3691 switch (fsType) {
64 caltinay 3750 case DegreesOfFreedom:
65     case ReducedDegreesOfFreedom:
66 caltinay 3670 case Nodes:
67     case ReducedNodes:
68     case Elements:
69     case ReducedElements:
70     case FaceElements:
71     case ReducedFaceElements:
72     case Points:
73 caltinay 3691 return true;
74     default:
75 caltinay 3670 break;
76     }
77 caltinay 3691 return false;
78 caltinay 3670 }
79    
80 caltinay 3691 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
81 caltinay 3670 {
82 caltinay 3691 switch (fsType) {
83 caltinay 3750 case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84     case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85 caltinay 3691 case Nodes: return "Ripley_Nodes";
86 caltinay 3764 case ReducedNodes: return "Ripley_ReducedNodes";
87 caltinay 3691 case Elements: return "Ripley_Elements";
88 caltinay 3764 case ReducedElements: return "Ripley_ReducedElements";
89     case FaceElements: return "Ripley_FaceElements";
90     case ReducedFaceElements: return "Ripley_ReducedFaceElements";
91 caltinay 3691 case Points: return "Ripley_Points";
92     default:
93 caltinay 3670 break;
94     }
95 caltinay 3691 return "Invalid function space type code";
96 caltinay 3670 }
97    
98 caltinay 3691 pair<int,int> RipleyDomain::getDataShape(int fsType) const
99 caltinay 3670 {
100 caltinay 3691 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
101     switch (fsType) {
102 caltinay 3670 case Nodes:
103 caltinay 3766 return pair<int,int>(1, getNumNodes());
104 caltinay 3744 case ReducedNodes: //FIXME: reduced
105 caltinay 3766 if (getCoarseMesh())
106     return getCoarseMesh()->getDataShape(Nodes);
107     break;
108 caltinay 3750 case DegreesOfFreedom:
109     case ReducedDegreesOfFreedom: //FIXME: reduced
110     return pair<int,int>(1, getNumDOF());
111 caltinay 3670 case Elements:
112 caltinay 3691 return pair<int,int>(ptsPerSample, getNumElements());
113 caltinay 3670 case FaceElements:
114 caltinay 3691 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
115 caltinay 3711 case ReducedElements:
116     return pair<int,int>(1, getNumElements());
117     case ReducedFaceElements:
118     return pair<int,int>(1, getNumFaceElements());
119 caltinay 3697 case Points:
120 caltinay 3733 return pair<int,int>(1, 0); //FIXME: dirac
121 caltinay 3670 default:
122     break;
123     }
124    
125 caltinay 3691 stringstream msg;
126 caltinay 3740 msg << "getDataShape(): Unsupported function space type " << fsType
127     << " for " << getDescription();
128 caltinay 3691 throw RipleyException(msg.str());
129 caltinay 3670 }
130    
131 caltinay 3691 string RipleyDomain::showTagNames() const
132 caltinay 3670 {
133 caltinay 3691 stringstream ret;
134     TagMap::const_iterator it;
135     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
136     if (it!=m_tagMap.begin()) ret << ", ";
137     ret << it->first;
138 caltinay 3670 }
139 caltinay 3691 return ret.str();
140 caltinay 3670 }
141    
142     bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
143     {
144 caltinay 3691 /*
145     The idea is to use equivalence classes (i.e. types which can be
146     interpolated back and forth):
147 caltinay 3750 class 0: DOF <-> Nodes
148     class 1: ReducedDOF <-> ReducedNodes
149 caltinay 3740 class 2: Points
150     class 3: Elements
151     class 4: ReducedElements
152     class 5: FaceElements
153     class 6: ReducedFaceElements
154 caltinay 3670
155 caltinay 3691 There is also a set of lines. Interpolation is possible down a line but not
156     between lines.
157 caltinay 3740 class 0 and 1 belong to all lines so aren't considered.
158     line 0: class 2
159     line 1: class 3,4
160     line 2: class 5,6
161 caltinay 3750
162     For classes with multiple members (eg class 1) we have vars to record if
163     there is at least one instance. e.g. hasnodes is true if we have at least
164     one instance of Nodes.
165 caltinay 3670 */
166     if (fs.empty())
167     return false;
168 caltinay 3740 vector<bool> hasclass(7, false);
169     vector<int> hasline(3, 0);
170 caltinay 3750 bool hasnodes=false;
171     bool hasrednodes=false;
172 caltinay 3740 for (size_t i=0; i<fs.size(); ++i) {
173 caltinay 3670 switch (fs[i]) {
174 caltinay 3750 case Nodes: hasnodes=true; // fall through
175     case DegreesOfFreedom:
176 caltinay 3740 hasclass[0]=true;
177 caltinay 3670 break;
178 caltinay 3750 case ReducedNodes: hasrednodes=true; // fall through
179     case ReducedDegreesOfFreedom:
180 caltinay 3740 hasclass[1]=true;
181 caltinay 3670 break;
182     case Points:
183     hasline[0]=1;
184 caltinay 3740 hasclass[2]=true;
185 caltinay 3670 break;
186     case Elements:
187 caltinay 3750 hasclass[3]=true;
188 caltinay 3670 hasline[1]=1;
189     break;
190     case ReducedElements:
191 caltinay 3740 hasclass[4]=true;
192 caltinay 3670 hasline[1]=1;
193     break;
194     case FaceElements:
195 caltinay 3750 hasclass[5]=true;
196 caltinay 3670 hasline[2]=1;
197     break;
198     case ReducedFaceElements:
199 caltinay 3740 hasclass[6]=true;
200 caltinay 3670 hasline[2]=1;
201     break;
202     default:
203     return false;
204     }
205     }
206 caltinay 3740 int numLines=hasline[0]+hasline[1]+hasline[2];
207 caltinay 3670
208     // fail if we have more than one leaf group
209     // = there are at least two branches we can't interpolate between
210 caltinay 3740 if (numLines > 1)
211 caltinay 3670 return false;
212 caltinay 3740 else if (numLines==1) {
213 caltinay 3750 // we have points
214 caltinay 3691 if (hasline[0]==1)
215 caltinay 3670 resultcode=Points;
216 caltinay 3691 else if (hasline[1]==1) {
217 caltinay 3740 if (hasclass[4])
218 caltinay 3670 resultcode=ReducedElements;
219 caltinay 3691 else
220 caltinay 3670 resultcode=Elements;
221 caltinay 3740 } else { // hasline[2]==1
222     if (hasclass[6])
223 caltinay 3670 resultcode=ReducedFaceElements;
224 caltinay 3691 else
225 caltinay 3670 resultcode=FaceElements;
226 caltinay 3740 }
227     } else { // numLines==0
228     if (hasclass[1])
229 caltinay 3750 // something from class 1
230     resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
231     else // something from class 0
232     resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
233 caltinay 3670 }
234     return true;
235     }
236    
237 caltinay 3732 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
238     int fsType_target) const
239     {
240 caltinay 3764 if (!isValidFunctionSpaceType(fsType_target)) {
241 caltinay 3732 stringstream msg;
242     msg << "probeInterpolationOnDomain(): Invalid functionspace type "
243 caltinay 3740 << fsType_target << " for " << getDescription();
244 caltinay 3732 throw RipleyException(msg.str());
245     }
246    
247     switch (fsType_source) {
248     case Nodes:
249 caltinay 3750 case DegreesOfFreedom:
250 caltinay 3732 return true;
251     case ReducedNodes:
252 caltinay 3750 case ReducedDegreesOfFreedom:
253     return (fsType_target != Nodes &&
254     fsType_target != DegreesOfFreedom);
255 caltinay 3732 case Elements:
256     return (fsType_target==Elements ||
257     fsType_target==ReducedElements);
258     case FaceElements:
259     return (fsType_target==FaceElements ||
260     fsType_target==ReducedFaceElements);
261 caltinay 3764 case ReducedElements:
262 caltinay 3732 case ReducedFaceElements:
263     case Points:
264 caltinay 3764 return (fsType_target==fsType_source);
265 caltinay 3732
266     default: {
267     stringstream msg;
268     msg << "probeInterpolationOnDomain(): Invalid functionspace type "
269 caltinay 3740 << fsType_source << " for " << getDescription();
270 caltinay 3732 throw RipleyException(msg.str());
271     }
272     }
273     }
274    
275 caltinay 3701 void RipleyDomain::interpolateOnDomain(escript::Data& target,
276     const escript::Data& in) const
277     {
278     const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
279     const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
280     if (inDomain != *this)
281     throw RipleyException("Illegal domain of interpolant");
282     if (targetDomain != *this)
283     throw RipleyException("Illegal domain of interpolation target");
284    
285     stringstream msg;
286     msg << "interpolateOnDomain() not implemented for function space "
287     << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
288     << " -> "
289     << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
290    
291 caltinay 3744 const int inFS = in.getFunctionSpace().getTypeCode();
292     const int outFS = target.getFunctionSpace().getTypeCode();
293 caltinay 3701
294 caltinay 3744 // simplest case: 1:1 copy
295     if (inFS==outFS) {
296     copyData(target, *const_cast<escript::Data*>(&in));
297     // not allowed: reduced->non-reduced
298 caltinay 3750 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
299     && (outFS==Nodes || outFS==DegreesOfFreedom)) {
300 caltinay 3744 throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
301     } else if ((inFS==Elements && outFS==ReducedElements)
302     || (inFS==FaceElements && outFS==ReducedFaceElements)) {
303     averageData(target, *const_cast<escript::Data*>(&in));
304     } else {
305     switch (inFS) {
306     case Nodes:
307 caltinay 3750 case ReducedNodes: //FIXME: reduced
308 caltinay 3744 switch (outFS) {
309 caltinay 3750 case DegreesOfFreedom:
310     case ReducedDegreesOfFreedom: //FIXME: reduced
311 caltinay 3755 if (getMPISize()==1)
312     copyData(target, *const_cast<escript::Data*>(&in));
313     else
314     nodesToDOF(target,*const_cast<escript::Data*>(&in));
315 caltinay 3744 break;
316 caltinay 3701
317 caltinay 3755 case Nodes:
318 caltinay 3766 copyData(target, *const_cast<escript::Data*>(&in));
319 caltinay 3755 case ReducedNodes: //FIXME: reduced
320     break;
321    
322 caltinay 3744 case Elements:
323     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
324     break;
325 caltinay 3711
326 caltinay 3744 case ReducedElements:
327     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
328     break;
329 caltinay 3701
330 caltinay 3744 case FaceElements:
331     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
332     break;
333 caltinay 3711
334 caltinay 3744 case ReducedFaceElements:
335     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
336     break;
337    
338     default:
339     throw RipleyException(msg.str());
340     }
341     break;
342 caltinay 3750
343     case DegreesOfFreedom:
344     case ReducedDegreesOfFreedom: //FIXME: reduced
345 caltinay 3755 switch (outFS) {
346     case Nodes:
347     case ReducedNodes: //FIXME: reduced
348     if (getMPISize()==1)
349     copyData(target, *const_cast<escript::Data*>(&in));
350     else
351     dofToNodes(target, *const_cast<escript::Data*>(&in));
352     break;
353    
354     case DegreesOfFreedom:
355     case ReducedDegreesOfFreedom: //FIXME: reduced
356     copyData(target, *const_cast<escript::Data*>(&in));
357     break;
358    
359     case Elements:
360     case ReducedElements:
361     if (getMPISize()==1) {
362     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
363     } else {
364 caltinay 3764 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
365     escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
366 caltinay 3755 interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
367     }
368     break;
369    
370     case FaceElements:
371     case ReducedFaceElements:
372     if (getMPISize()==1) {
373     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
374     } else {
375 caltinay 3764 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
376     escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
377     interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
378 caltinay 3755 }
379     break;
380    
381     default:
382     throw RipleyException(msg.str());
383     }
384 caltinay 3754 break;
385 caltinay 3744 default:
386     throw RipleyException(msg.str());
387     }
388 caltinay 3701 }
389     }
390    
391 caltinay 3691 escript::Data RipleyDomain::getX() const
392 caltinay 3670 {
393 caltinay 3691 return escript::continuousFunction(*this).getX();
394     }
395    
396     escript::Data RipleyDomain::getNormal() const
397     {
398     return escript::functionOnBoundary(*this).getNormal();
399     }
400    
401     escript::Data RipleyDomain::getSize() const
402     {
403     return escript::function(*this).getSize();
404     }
405    
406     void RipleyDomain::setToX(escript::Data& arg) const
407     {
408     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
409     *(arg.getFunctionSpace().getDomain()));
410     if (argDomain != *this)
411     throw RipleyException("setToX: Illegal domain of data point locations");
412     if (!arg.isExpanded())
413     throw RipleyException("setToX: Expanded Data object expected");
414    
415     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
416     assembleCoordinates(arg);
417     } else {
418     // interpolate the result
419     escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
420     assembleCoordinates(contData);
421     interpolateOnDomain(arg, contData);
422 caltinay 3670 }
423 caltinay 3691 }
424 caltinay 3670
425 caltinay 3764 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
426     {
427     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
428     *(arg.getFunctionSpace().getDomain()));
429     if (argDomain != *this)
430     throw RipleyException("setToGradient: Illegal domain of gradient argument");
431     const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
432     *(grad.getFunctionSpace().getDomain()));
433     if (gradDomain != *this)
434     throw RipleyException("setToGradient: Illegal domain of gradient");
435    
436     switch (grad.getFunctionSpace().getTypeCode()) {
437     case Elements:
438     case ReducedElements:
439     case FaceElements:
440     case ReducedFaceElements:
441     break;
442     default: {
443     stringstream msg;
444     msg << "setToGradient(): not supported for "
445     << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
446     throw RipleyException(msg.str());
447     }
448     }
449    
450     if (getMPISize()>1) {
451     if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
452     escript::Data contArg(arg, escript::continuousFunction(*this));
453     assembleGradient(grad, contArg);
454     } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
455     escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
456     assembleGradient(grad, contArg);
457     } else {
458     assembleGradient(grad, *const_cast<escript::Data*>(&arg));
459     }
460     } else {
461     assembleGradient(grad, *const_cast<escript::Data*>(&arg));
462     }
463     }
464    
465     void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
466     {
467     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
468     *(arg.getFunctionSpace().getDomain()));
469     if (argDomain != *this)
470     throw RipleyException("setToIntegrals: Illegal domain of integration kernel");
471    
472     switch (arg.getFunctionSpace().getTypeCode()) {
473     case Nodes:
474     case ReducedNodes:
475     case DegreesOfFreedom:
476     case ReducedDegreesOfFreedom:
477     {
478     escript::Data funcArg(arg, escript::function(*this));
479     assembleIntegrate(integrals, funcArg);
480     }
481     break;
482     case Elements:
483     case ReducedElements:
484     case FaceElements:
485     case ReducedFaceElements:
486     assembleIntegrate(integrals, *const_cast<escript::Data*>(&arg));
487     break;
488     default: {
489     stringstream msg;
490     msg << "setToIntegrals(): not supported for "
491     << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
492     throw RipleyException(msg.str());
493     }
494     }
495    
496     }
497    
498 caltinay 3691 bool RipleyDomain::isCellOriented(int fsType) const
499     {
500     switch(fsType) {
501 caltinay 3670 case Nodes:
502 caltinay 3764 case ReducedNodes:
503 caltinay 3750 case DegreesOfFreedom:
504     case ReducedDegreesOfFreedom:
505 caltinay 3691 return false;
506 caltinay 3670 case Elements:
507 caltinay 3691 case FaceElements:
508     case Points:
509 caltinay 3670 case ReducedElements:
510     case ReducedFaceElements:
511 caltinay 3691 return true;
512 caltinay 3670 default:
513 caltinay 3691 break;
514 caltinay 3670 }
515 caltinay 3691 stringstream msg;
516 caltinay 3740 msg << "isCellOriented(): Illegal function space type " << fsType
517     << " on " << getDescription();
518 caltinay 3691 throw RipleyException(msg.str());
519 caltinay 3670 }
520    
521 caltinay 3722 bool RipleyDomain::canTag(int fsType) const
522     {
523     switch(fsType) {
524     case Nodes:
525     case Elements:
526     case ReducedElements:
527     case FaceElements:
528     case ReducedFaceElements:
529     return true;
530 caltinay 3750 case DegreesOfFreedom:
531     case ReducedDegreesOfFreedom:
532 caltinay 3722 case Points:
533 caltinay 3764 case ReducedNodes:
534 caltinay 3722 return false;
535     default:
536     break;
537     }
538     stringstream msg;
539 caltinay 3740 msg << "canTag(): Illegal function space type " << fsType << " on "
540     << getDescription();
541 caltinay 3722 throw RipleyException(msg.str());
542     }
543    
544     void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
545     {
546     IndexVector* target=NULL;
547     dim_t num=0;
548    
549     switch(fsType) {
550     case Nodes:
551     num=getNumNodes();
552     target=&m_nodeTags;
553     break;
554     case Elements:
555     case ReducedElements:
556     num=getNumElements();
557     target=&m_elementTags;
558     break;
559     case FaceElements:
560     case ReducedFaceElements:
561     num=getNumFaceElements();
562     target=&m_faceTags;
563     break;
564     default: {
565     stringstream msg;
566 caltinay 3740 msg << "setTags(): not implemented for "
567     << functionSpaceTypeAsString(fsType);
568 caltinay 3722 throw RipleyException(msg.str());
569     }
570     }
571     if (target->size() != num) {
572     target->assign(num, -1);
573     }
574    
575     escript::Data& mask=*const_cast<escript::Data*>(&cMask);
576     #pragma omp parallel for
577     for (index_t i=0; i<num; i++) {
578     if (mask.getSampleDataRO(i)[0] > 0) {
579     (*target)[i]=newTag;
580     }
581     }
582     updateTagsInUse(fsType);
583     }
584    
585     int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
586     {
587     switch(fsType) {
588     case Nodes:
589     if (m_nodeTags.size() > sampleNo)
590     return m_nodeTags[sampleNo];
591     break;
592     case Elements:
593     case ReducedElements:
594     if (m_elementTags.size() > sampleNo)
595     return m_elementTags[sampleNo];
596     break;
597     case FaceElements:
598     case ReducedFaceElements:
599     if (m_faceTags.size() > sampleNo)
600     return m_faceTags[sampleNo];
601     break;
602     default: {
603     stringstream msg;
604 caltinay 3740 msg << "getTagFromSampleNo(): not implemented for "
605     << functionSpaceTypeAsString(fsType);
606 caltinay 3722 throw RipleyException(msg.str());
607     }
608     }
609     return -1;
610     }
611    
612     int RipleyDomain::getNumberOfTagsInUse(int fsType) const
613     {
614     switch(fsType) {
615     case Nodes:
616     return m_nodeTagsInUse.size();
617     case Elements:
618     case ReducedElements:
619     return m_elementTagsInUse.size();
620     case FaceElements:
621     case ReducedFaceElements:
622     return m_faceTagsInUse.size();
623     default: {
624     stringstream msg;
625 caltinay 3740 msg << "getNumberOfTagsInUse(): not implemented for "
626     << functionSpaceTypeAsString(fsType);
627 caltinay 3722 throw RipleyException(msg.str());
628     }
629     }
630     }
631    
632     const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
633     {
634     switch(fsType) {
635     case Nodes:
636     return &m_nodeTagsInUse[0];
637     case Elements:
638     case ReducedElements:
639     return &m_elementTagsInUse[0];
640     case FaceElements:
641     case ReducedFaceElements:
642     return &m_faceTagsInUse[0];
643     default: {
644     stringstream msg;
645 caltinay 3740 msg << "borrowListOfTagsInUse(): not implemented for "
646     << functionSpaceTypeAsString(fsType);
647 caltinay 3722 throw RipleyException(msg.str());
648     }
649     }
650     }
651    
652 caltinay 3691 void RipleyDomain::Print_Mesh_Info(const bool full) const
653 caltinay 3670 {
654 caltinay 3691 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
655     << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
656     cout << "Number of dimensions: " << m_numDim << endl;
657 caltinay 3670
658 caltinay 3691 // write tags
659     if (m_tagMap.size() > 0) {
660     cout << "Tags:" << endl;
661     TagMap::const_iterator it;
662     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
663     cout << " " << setw(5) << it->second << " "
664     << it->first << endl;
665     }
666 caltinay 3670 }
667     }
668    
669     int RipleyDomain::getSystemMatrixTypeId(const int solver,
670     const int preconditioner, const int package, const bool symmetry) const
671     {
672 caltinay 3691 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
673     package, symmetry, m_mpiInfo);
674 caltinay 3670 }
675    
676     int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
677     const int package, const bool symmetry) const
678     {
679 caltinay 3691 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
680     package, symmetry, m_mpiInfo);
681 caltinay 3670 }
682    
683 caltinay 3691 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
684     const escript::FunctionSpace& row_functionspace,
685     const int column_blocksize,
686     const escript::FunctionSpace& column_functionspace,
687     const int type) const
688 caltinay 3670 {
689 caltinay 3691 bool reduceRowOrder=false;
690     bool reduceColOrder=false;
691     // is the domain right?
692     const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
693 caltinay 3764 if (row_domain != *this)
694 caltinay 3740 throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
695 caltinay 3691 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
696 caltinay 3764 if (col_domain != *this)
697 caltinay 3740 throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
698 caltinay 3691 // is the function space type right?
699 caltinay 3750 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
700 caltinay 3691 reduceRowOrder=true;
701 caltinay 3750 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
702 caltinay 3740 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
703 caltinay 3750 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
704 caltinay 3691 reduceColOrder=true;
705 caltinay 3750 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
706 caltinay 3740 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
707 caltinay 3691
708     // generate matrix
709     Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
710     Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
711     row_blocksize, column_blocksize, FALSE);
712     paso::checkPasoError();
713     Paso_SystemMatrixPattern_free(pattern);
714     escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
715     row_functionspace, column_blocksize, column_functionspace));
716     return sma;
717 caltinay 3670 }
718    
719 caltinay 3748 void RipleyDomain::addPDEToSystem(
720     escript::AbstractSystemMatrix& mat, escript::Data& rhs,
721     const escript::Data& A, const escript::Data& B, const escript::Data& C,
722     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
723     const escript::Data& d, const escript::Data& y,
724     const escript::Data& d_contact, const escript::Data& y_contact,
725     const escript::Data& d_dirac,const escript::Data& y_dirac) const
726     {
727     if (!d_contact.isEmpty() || !y_contact.isEmpty())
728     throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
729    
730     paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
731     if (!sma)
732     throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
733    
734     if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
735     throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
736    
737 caltinay 3761 int fsType=UNKNOWN;
738     fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
739     fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
740     fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
741     fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
742     fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
743     fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
744     fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
745     fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
746     fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
747     fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
748     if (fsType==UNKNOWN)
749     return;
750 caltinay 3748
751 caltinay 3761 const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
752    
753     //TODO: more input param checks (shape, etc)
754    
755 caltinay 3748 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
756    
757     if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
758     throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
759    
760     const int numEq=S->logical_row_block_size;
761     const int numComp=S->logical_col_block_size;
762    
763     if (numEq != numComp)
764     throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
765     //TODO: more system matrix checks
766    
767     if (numEq==1)
768 caltinay 3761 if (reducedOrder)
769     assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y, d, y);
770     else
771     assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
772 caltinay 3748 else
773 caltinay 3761 if (reducedOrder)
774     assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y, d, y);
775     else
776     assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
777 caltinay 3748 }
778    
779 caltinay 3760 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
780     const escript::Data& Y, const escript::Data& y,
781     const escript::Data& y_contact, const escript::Data& y_dirac) const
782     {
783     if (!y_contact.isEmpty())
784     throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
785    
786     if (rhs.isEmpty()) {
787     if (!X.isEmpty() || !Y.isEmpty())
788     throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
789     else
790     return;
791     }
792    
793     if (rhs.getDataPointSize() == 1)
794     assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
795     else
796     assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
797     }
798    
799 caltinay 3697 void RipleyDomain::setNewX(const escript::Data& arg)
800     {
801     throw RipleyException("setNewX(): Operation not supported");
802     }
803    
804 caltinay 3722 //protected
805 caltinay 3744 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
806 caltinay 3722 {
807     const dim_t numComp = in.getDataPointSize();
808 caltinay 3744 const dim_t dpp = in.getNumDataPointsPerSample();
809 caltinay 3740 out.requireWrite();
810 caltinay 3722 #pragma omp parallel for
811     for (index_t i=0; i<in.getNumSamples(); i++) {
812     const double* src = in.getSampleDataRO(i);
813 caltinay 3744 double* dest = out.getSampleDataRW(i);
814     for (index_t c=0; c<numComp; c++) {
815     double res=0.;
816     for (index_t q=0; q<dpp; q++)
817     res+=src[c+q*numComp];
818     *dest++ = res/dpp;
819     }
820     }
821     }
822    
823     //protected
824     void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
825     {
826     const dim_t numComp = in.getDataPointSize();
827     out.requireWrite();
828     #pragma omp parallel for
829     for (index_t i=0; i<in.getNumSamples(); i++) {
830     const double* src = in.getSampleDataRO(i);
831 caltinay 3722 copy(src, src+numComp, out.getSampleDataRW(i));
832     }
833     }
834    
835     //protected
836     void RipleyDomain::updateTagsInUse(int fsType) const
837     {
838     IndexVector* tagsInUse=NULL;
839     const IndexVector* tags=NULL;
840     switch(fsType) {
841     case Nodes:
842     tags=&m_nodeTags;
843     tagsInUse=&m_nodeTagsInUse;
844     break;
845     case Elements:
846     case ReducedElements:
847     tags=&m_elementTags;
848     tagsInUse=&m_elementTagsInUse;
849     break;
850     case FaceElements:
851     case ReducedFaceElements:
852     tags=&m_faceTags;
853     tagsInUse=&m_faceTagsInUse;
854     break;
855     default:
856     return;
857     }
858    
859     // gather global unique tag values from tags into tagsInUse
860     tagsInUse->clear();
861     index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
862    
863     while (true) {
864     // find smallest value bigger than lastFoundValue
865     minFoundValue = INDEX_T_MAX;
866     #pragma omp parallel private(local_minFoundValue)
867     {
868     local_minFoundValue = minFoundValue;
869     #pragma omp for schedule(static) nowait
870     for (size_t i = 0; i < tags->size(); i++) {
871     const index_t v = (*tags)[i];
872     if ((v > lastFoundValue) && (v < local_minFoundValue))
873     local_minFoundValue = v;
874     }
875     #pragma omp critical
876     {
877     if (local_minFoundValue < minFoundValue)
878     minFoundValue = local_minFoundValue;
879     }
880     }
881     #ifdef ESYS_MPI
882     local_minFoundValue = minFoundValue;
883     MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
884     #endif
885 caltinay 3740
886 caltinay 3722 // if we found a new value add it to the tagsInUse vector
887     if (minFoundValue < INDEX_T_MAX) {
888     tagsInUse->push_back(minFoundValue);
889     lastFoundValue = minFoundValue;
890     } else
891     break;
892     }
893     }
894    
895 caltinay 3756 //protected
896     Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
897     const IndexVector& index, const dim_t M, const dim_t N) const
898     {
899     // paso will manage the memory
900     index_t* indexC = MEMALLOC(index.size(), index_t);
901     index_t* ptrC = MEMALLOC(ptr.size(), index_t);
902     copy(index.begin(), index.end(), indexC);
903     copy(ptr.begin(), ptr.end(), ptrC);
904     return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
905     }
906    
907     //protected
908     Paso_Pattern* RipleyDomain::createMainPattern() const
909     {
910     IndexVector ptr(1,0);
911     IndexVector index;
912    
913     for (index_t i=0; i<getNumDOF(); i++) {
914     // add the DOF itself
915     index.push_back(i);
916     const dim_t num=insertNeighbourNodes(index, i);
917     ptr.push_back(ptr.back()+num+1);
918     }
919    
920     return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
921     }
922    
923     //protected
924     void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
925     const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
926     {
927     IndexVector ptr(1,0);
928     IndexVector index;
929     for (index_t i=0; i<getNumDOF(); i++) {
930     index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
931     ptr.push_back(ptr.back()+colIndices[i].size());
932     }
933    
934     const dim_t M=ptr.size()-1;
935     *colPattern=createPasoPattern(ptr, index, M, N);
936    
937     IndexVector rowPtr(1,0);
938     IndexVector rowIndex;
939     for (dim_t id=0; id<N; id++) {
940     dim_t n=0;
941     for (dim_t i=0; i<M; i++) {
942     for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
943     if (index[j]==id) {
944     rowIndex.push_back(i);
945     n++;
946     break;
947     }
948     }
949     }
950     rowPtr.push_back(rowPtr.back()+n);
951     }
952    
953     // M and N are now reversed
954     *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
955     }
956    
957     //protected
958     void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
959     const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
960     dim_t num_Sol, const vector<double>& array) const
961     {
962     const dim_t numMyCols = mat->pattern->mainPattern->numInput;
963     const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
964 caltinay 3760 const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
965     const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
966 caltinay 3756
967     const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
968     const index_t* mainBlock_index = mat->mainBlock->pattern->index;
969     double* mainBlock_val = mat->mainBlock->val;
970     const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
971     const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
972     double* col_coupleBlock_val = mat->col_coupleBlock->val;
973     const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
974     const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
975     double* row_coupleBlock_val = mat->row_coupleBlock->val;
976    
977     for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
978     // down columns of array
979 caltinay 3760 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
980     const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
981     // only look at the matrix rows stored on this processor
982     if (i_row < numMyRows) {
983     for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
984     for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
985     const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
986     if (i_col < numMyCols) {
987     for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
988     if (mainBlock_index[k] == i_col) {
989     for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
990     const dim_t i_Sol=ic+mat->col_block_size*l_col;
991     for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
992     const dim_t i_Eq=ir+mat->row_block_size*l_row;
993     mainBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
994     }
995     }
996     break;
997     }
998     }
999     } else {
1000     for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1001     if (col_coupleBlock_index[k] == i_col - numMyCols) {
1002     for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1003     const dim_t i_Sol=ic+mat->col_block_size*l_col;
1004     for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1005     const dim_t i_Eq=ir+mat->row_block_size*l_row;
1006     col_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1007     }
1008     }
1009     break;
1010     }
1011     }
1012 caltinay 3756 }
1013     }
1014     }
1015 caltinay 3760 } else {
1016     for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1017     // across rows of array
1018     for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1019     const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1020     if (i_col < numMyCols) {
1021     for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1022     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1023     {
1024     if (row_coupleBlock_index[k] == i_col) {
1025     for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1026     const dim_t i_Sol=ic+mat->col_block_size*l_col;
1027     for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1028     const dim_t i_Eq=ir+mat->row_block_size*l_row;
1029     row_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1030     }
1031     }
1032     break;
1033     }
1034     }
1035 caltinay 3756 }
1036     }
1037     }
1038     }
1039     }
1040     }
1041     }
1042    
1043 caltinay 3691 //
1044     // the methods that follow have to be implemented by the subclasses
1045     //
1046    
1047     string RipleyDomain::getDescription() const
1048 caltinay 3670 {
1049 caltinay 3691 throw RipleyException("getDescription() not implemented");
1050 caltinay 3670 }
1051    
1052 caltinay 3691 void RipleyDomain::write(const string& filename) const
1053 caltinay 3670 {
1054 caltinay 3691 throw RipleyException("write() not implemented");
1055 caltinay 3670 }
1056    
1057 caltinay 3691 void RipleyDomain::dump(const string& filename) const
1058 caltinay 3670 {
1059 caltinay 3691 throw RipleyException("dump() not implemented");
1060 caltinay 3670 }
1061    
1062 caltinay 3691 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
1063     {
1064     throw RipleyException("borrowSampleReferenceIDs() not implemented");
1065     }
1066 caltinay 3670
1067 caltinay 3691 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1068 caltinay 3670 {
1069 caltinay 3691 throw RipleyException("interpolateACross() not implemented");
1070 caltinay 3670 }
1071    
1072 caltinay 3691 bool RipleyDomain::probeInterpolationACross(int fsType_source,
1073     const escript::AbstractDomain&, int fsType_target) const
1074 caltinay 3670 {
1075 caltinay 3691 throw RipleyException("probeInterpolationACross() not implemented");
1076 caltinay 3670 }
1077    
1078 caltinay 3691 void RipleyDomain::setToNormal(escript::Data& normal) const
1079 caltinay 3670 {
1080 caltinay 3691 throw RipleyException("setToNormal() not implemented");
1081 caltinay 3670 }
1082    
1083 caltinay 3691 void RipleyDomain::setToSize(escript::Data& size) const
1084 caltinay 3670 {
1085 caltinay 3691 throw RipleyException("setToSize() not implemented");
1086 caltinay 3670 }
1087    
1088 caltinay 3691 bool RipleyDomain::ownSample(int fsType, index_t id) const
1089 caltinay 3670 {
1090 caltinay 3691 throw RipleyException("ownSample() not implemented");
1091 caltinay 3670 }
1092    
1093 caltinay 3691 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1094     const escript::Data& D, const escript::Data& d,
1095     const escript::Data& d_dirac, const bool useHRZ) const
1096     {
1097     throw RipleyException("addPDEToLumpedSystem() not implemented");
1098 caltinay 3670 }
1099    
1100 caltinay 3691 void RipleyDomain::addPDEToTransportProblem(
1101     escript::AbstractTransportProblem& tp,
1102     escript::Data& source, const escript::Data& M,
1103     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1104     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1105     const escript::Data& d, const escript::Data& y,
1106     const escript::Data& d_contact, const escript::Data& y_contact,
1107     const escript::Data& d_dirac, const escript::Data& y_dirac) const
1108 caltinay 3670 {
1109 caltinay 3691 throw RipleyException("addPDEToTransportProblem() not implemented");
1110 caltinay 3670 }
1111    
1112 caltinay 3691 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
1113     const int blocksize, const escript::FunctionSpace& functionspace,
1114     const int type) const
1115 caltinay 3670 {
1116 caltinay 3691 throw RipleyException("newTransportProblem() not implemented");
1117 caltinay 3670 }
1118    
1119 caltinay 3691 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1120     bool reducedColOrder) const
1121 caltinay 3670 {
1122 caltinay 3691 throw RipleyException("getPattern() not implemented");
1123 caltinay 3670 }
1124    
1125 caltinay 3691 dim_t RipleyDomain::getNumDataPointsGlobal() const
1126 caltinay 3670 {
1127 caltinay 3691 throw RipleyException("getNumDataPointsGlobal() not implemented");
1128 caltinay 3670 }
1129    
1130 caltinay 3697 IndexVector RipleyDomain::getNumNodesPerDim() const
1131     {
1132     throw RipleyException("getNumNodesPerDim() not implemented");
1133     }
1134    
1135     IndexVector RipleyDomain::getNumElementsPerDim() const
1136     {
1137     throw RipleyException("getNumElementsPerDim() not implemented");
1138     }
1139    
1140     IndexVector RipleyDomain::getNumFacesPerBoundary() const
1141     {
1142     throw RipleyException("getNumFacesPerBoundary() not implemented");
1143     }
1144    
1145     IndexVector RipleyDomain::getNodeDistribution() const
1146     {
1147     throw RipleyException("getNodeDistribution() not implemented");
1148     }
1149    
1150 caltinay 3766 IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1151     {
1152     throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1153     }
1154    
1155 caltinay 3697 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1156     {
1157     throw RipleyException("getFirstCoordAndSpacing() not implemented");
1158     }
1159    
1160 caltinay 3691 dim_t RipleyDomain::getNumFaceElements() const
1161 caltinay 3670 {
1162 caltinay 3691 throw RipleyException("getNumFaceElements() not implemented");
1163 caltinay 3670 }
1164    
1165 caltinay 3691 dim_t RipleyDomain::getNumElements() const
1166 caltinay 3670 {
1167 caltinay 3691 throw RipleyException("getNumElements() not implemented");
1168 caltinay 3670 }
1169    
1170 caltinay 3691 dim_t RipleyDomain::getNumNodes() const
1171 caltinay 3670 {
1172 caltinay 3691 throw RipleyException("getNumNodes() not implemented");
1173 caltinay 3670 }
1174    
1175 caltinay 3750 dim_t RipleyDomain::getNumDOF() const
1176     {
1177     throw RipleyException("getNumDOF() not implemented");
1178     }
1179    
1180 caltinay 3766 escript::Domain_ptr RipleyDomain::getCoarseMesh() const
1181     {
1182     throw RipleyException("getCoarseMesh() not implemented");
1183     }
1184    
1185 caltinay 3756 dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1186     {
1187     throw RipleyException("insertNeighbourNodes() not implemented");
1188     }
1189    
1190 caltinay 3691 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1191 caltinay 3670 {
1192 caltinay 3691 throw RipleyException("assembleCoordinates() not implemented");
1193 caltinay 3670 }
1194    
1195 caltinay 3764 void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1196     {
1197     throw RipleyException("assembleGradient() not implemented");
1198     }
1199    
1200     void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1201     {
1202     throw RipleyException("assembleIntegrate() not implemented");
1203     }
1204    
1205 caltinay 3748 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1206     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1207     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1208     const escript::Data& d, const escript::Data& y) const
1209     {
1210     throw RipleyException("assemblePDESingle() not implemented");
1211     }
1212    
1213 caltinay 3761 void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1214     escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1215     const escript::Data& C, const escript::Data& D, const escript::Data& X,
1216     const escript::Data& Y, const escript::Data& d, const escript::Data& y) const
1217     {
1218     throw RipleyException("assemblePDESingleReduced() not implemented");
1219     }
1220    
1221 caltinay 3748 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1222     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1223     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1224     const escript::Data& d, const escript::Data& y) const
1225     {
1226     throw RipleyException("assemblePDESystem() not implemented");
1227     }
1228    
1229 caltinay 3761 void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1230     escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1231     const escript::Data& C, const escript::Data& D, const escript::Data& X,
1232     const escript::Data& Y, const escript::Data& d, const escript::Data& y) const
1233     {
1234     throw RipleyException("assemblePDESystemReduced() not implemented");
1235     }
1236    
1237 caltinay 3711 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1238 caltinay 3701 {
1239     throw RipleyException("interpolateNodesOnElements() not implemented");
1240     }
1241 caltinay 3691
1242 caltinay 3711 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1243 caltinay 3701 {
1244     throw RipleyException("interpolateNodesOnFaces() not implemented");
1245     }
1246    
1247 caltinay 3750 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1248     {
1249     throw RipleyException("nodesToDOF() not implemented");
1250     }
1251 caltinay 3701
1252 caltinay 3754 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1253     {
1254     throw RipleyException("dofToNodes() not implemented");
1255     }
1256    
1257 caltinay 3670 } // end of namespace ripley
1258    

  ViewVC Help
Powered by ViewVC 1.1.26