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

Annotation of /branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3776 - (hide annotations)
Thu Jan 19 03:48:35 2012 UTC (7 years, 9 months ago) by caltinay
File size: 45966 byte(s)
Finished 2D - all 2D PDE (non-lumped) tests pass now

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

  ViewVC Help
Powered by ViewVC 1.1.26