/[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 3785 - (hide annotations)
Wed Jan 25 04:00:33 2012 UTC (7 years, 9 months ago) by caltinay
File size: 45781 byte(s)
Fixed early deallocation of matrix pattern.

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 caltinay 3777 if (in.actsExpanded())
301     averageData(target, *const_cast<escript::Data*>(&in));
302     else
303     copyData(target, *const_cast<escript::Data*>(&in));
304 caltinay 3744 } 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 3776 case ReducedNodes: //FIXME: reduced
319 caltinay 3766 copyData(target, *const_cast<escript::Data*>(&in));
320 caltinay 3755 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     escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
714     row_functionspace, column_blocksize, column_functionspace));
715     return sma;
716 caltinay 3670 }
717    
718 caltinay 3748 void RipleyDomain::addPDEToSystem(
719     escript::AbstractSystemMatrix& mat, escript::Data& rhs,
720     const escript::Data& A, const escript::Data& B, const escript::Data& C,
721     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
722     const escript::Data& d, const escript::Data& y,
723     const escript::Data& d_contact, const escript::Data& y_contact,
724     const escript::Data& d_dirac,const escript::Data& y_dirac) const
725     {
726     if (!d_contact.isEmpty() || !y_contact.isEmpty())
727     throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
728    
729     paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
730     if (!sma)
731     throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
732    
733     if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
734     throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
735    
736 caltinay 3761 int fsType=UNKNOWN;
737     fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
738     fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
739     fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
740     fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
741     fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
742     fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
743     fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
744     fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
745     fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
746     fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
747     if (fsType==UNKNOWN)
748     return;
749 caltinay 3748
750 caltinay 3761 const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
751    
752     //TODO: more input param checks (shape, etc)
753    
754 caltinay 3748 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
755    
756     if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
757     throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
758    
759     const int numEq=S->logical_row_block_size;
760     const int numComp=S->logical_col_block_size;
761    
762     if (numEq != numComp)
763     throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
764     //TODO: more system matrix checks
765    
766 caltinay 3769 if (numEq==1) {
767     if (reducedOrder) {
768     assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y);
769 caltinay 3776 assemblePDEBoundarySingleReduced(S, rhs, d, y);
770 caltinay 3769 } else {
771     assemblePDESingle(S, rhs, A, B, C, D, X, Y);
772 caltinay 3776 assemblePDEBoundarySingle(S, rhs, d, y);
773 caltinay 3769 }
774     } else {
775     if (reducedOrder) {
776     assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y);
777 caltinay 3776 assemblePDEBoundarySystemReduced(S, rhs, d, y);
778 caltinay 3769 } else {
779     assemblePDESystem(S, rhs, A, B, C, D, X, Y);
780 caltinay 3776 assemblePDEBoundarySystem(S, rhs, d, y);
781 caltinay 3769 }
782     }
783 caltinay 3748 }
784    
785 caltinay 3760 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
786     const escript::Data& Y, const escript::Data& y,
787     const escript::Data& y_contact, const escript::Data& y_dirac) const
788     {
789     if (!y_contact.isEmpty())
790     throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
791    
792     if (rhs.isEmpty()) {
793     if (!X.isEmpty() || !Y.isEmpty())
794     throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
795     else
796     return;
797     }
798    
799 caltinay 3769 if (rhs.getDataPointSize() == 1) {
800     assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
801     assemblePDEBoundarySingle(NULL, rhs, escript::Data(), y);
802     } else {
803     assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
804     assemblePDEBoundarySystem(NULL, rhs, escript::Data(), y);
805     }
806 caltinay 3760 }
807    
808 caltinay 3697 void RipleyDomain::setNewX(const escript::Data& arg)
809     {
810     throw RipleyException("setNewX(): Operation not supported");
811     }
812    
813 caltinay 3722 //protected
814 caltinay 3744 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
815 caltinay 3722 {
816     const dim_t numComp = in.getDataPointSize();
817 caltinay 3744 const dim_t dpp = in.getNumDataPointsPerSample();
818 caltinay 3740 out.requireWrite();
819 caltinay 3722 #pragma omp parallel for
820     for (index_t i=0; i<in.getNumSamples(); i++) {
821     const double* src = in.getSampleDataRO(i);
822 caltinay 3744 double* dest = out.getSampleDataRW(i);
823     for (index_t c=0; c<numComp; c++) {
824     double res=0.;
825     for (index_t q=0; q<dpp; q++)
826     res+=src[c+q*numComp];
827     *dest++ = res/dpp;
828     }
829     }
830     }
831    
832     //protected
833     void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
834     {
835     const dim_t numComp = in.getDataPointSize();
836     out.requireWrite();
837     #pragma omp parallel for
838     for (index_t i=0; i<in.getNumSamples(); i++) {
839     const double* src = in.getSampleDataRO(i);
840 caltinay 3722 copy(src, src+numComp, out.getSampleDataRW(i));
841     }
842     }
843    
844     //protected
845     void RipleyDomain::updateTagsInUse(int fsType) const
846     {
847     IndexVector* tagsInUse=NULL;
848     const IndexVector* tags=NULL;
849     switch(fsType) {
850     case Nodes:
851     tags=&m_nodeTags;
852     tagsInUse=&m_nodeTagsInUse;
853     break;
854     case Elements:
855     case ReducedElements:
856     tags=&m_elementTags;
857     tagsInUse=&m_elementTagsInUse;
858     break;
859     case FaceElements:
860     case ReducedFaceElements:
861     tags=&m_faceTags;
862     tagsInUse=&m_faceTagsInUse;
863     break;
864     default:
865     return;
866     }
867    
868     // gather global unique tag values from tags into tagsInUse
869     tagsInUse->clear();
870     index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
871    
872     while (true) {
873     // find smallest value bigger than lastFoundValue
874     minFoundValue = INDEX_T_MAX;
875     #pragma omp parallel private(local_minFoundValue)
876     {
877     local_minFoundValue = minFoundValue;
878     #pragma omp for schedule(static) nowait
879     for (size_t i = 0; i < tags->size(); i++) {
880     const index_t v = (*tags)[i];
881     if ((v > lastFoundValue) && (v < local_minFoundValue))
882     local_minFoundValue = v;
883     }
884     #pragma omp critical
885     {
886     if (local_minFoundValue < minFoundValue)
887     minFoundValue = local_minFoundValue;
888     }
889     }
890     #ifdef ESYS_MPI
891     local_minFoundValue = minFoundValue;
892     MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
893     #endif
894 caltinay 3740
895 caltinay 3722 // if we found a new value add it to the tagsInUse vector
896     if (minFoundValue < INDEX_T_MAX) {
897     tagsInUse->push_back(minFoundValue);
898     lastFoundValue = minFoundValue;
899     } else
900     break;
901     }
902     }
903    
904 caltinay 3756 //protected
905     Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
906     const IndexVector& index, const dim_t M, const dim_t N) const
907     {
908     // paso will manage the memory
909     index_t* indexC = MEMALLOC(index.size(), index_t);
910     index_t* ptrC = MEMALLOC(ptr.size(), index_t);
911     copy(index.begin(), index.end(), indexC);
912     copy(ptr.begin(), ptr.end(), ptrC);
913     return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
914     }
915    
916     //protected
917     Paso_Pattern* RipleyDomain::createMainPattern() const
918     {
919     IndexVector ptr(1,0);
920     IndexVector index;
921    
922     for (index_t i=0; i<getNumDOF(); i++) {
923     // add the DOF itself
924     index.push_back(i);
925     const dim_t num=insertNeighbourNodes(index, i);
926     ptr.push_back(ptr.back()+num+1);
927     }
928    
929     return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
930     }
931    
932     //protected
933     void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
934     const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
935     {
936     IndexVector ptr(1,0);
937     IndexVector index;
938     for (index_t i=0; i<getNumDOF(); i++) {
939     index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
940     ptr.push_back(ptr.back()+colIndices[i].size());
941     }
942    
943     const dim_t M=ptr.size()-1;
944     *colPattern=createPasoPattern(ptr, index, M, N);
945    
946     IndexVector rowPtr(1,0);
947     IndexVector rowIndex;
948     for (dim_t id=0; id<N; id++) {
949     dim_t n=0;
950     for (dim_t i=0; i<M; i++) {
951     for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
952     if (index[j]==id) {
953     rowIndex.push_back(i);
954     n++;
955     break;
956     }
957     }
958     }
959     rowPtr.push_back(rowPtr.back()+n);
960     }
961    
962     // M and N are now reversed
963     *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
964     }
965    
966     //protected
967     void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
968     const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
969     dim_t num_Sol, const vector<double>& array) const
970     {
971     const dim_t numMyCols = mat->pattern->mainPattern->numInput;
972     const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
973 caltinay 3760 const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
974     const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
975 caltinay 3756
976     const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
977     const index_t* mainBlock_index = mat->mainBlock->pattern->index;
978     double* mainBlock_val = mat->mainBlock->val;
979     const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
980     const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
981     double* col_coupleBlock_val = mat->col_coupleBlock->val;
982     const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
983     const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
984     double* row_coupleBlock_val = mat->row_coupleBlock->val;
985    
986     for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
987     // down columns of array
988 caltinay 3760 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
989     const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
990     // only look at the matrix rows stored on this processor
991     if (i_row < numMyRows) {
992     for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
993     for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
994     const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
995     if (i_col < numMyCols) {
996     for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
997     if (mainBlock_index[k] == i_col) {
998     for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
999     const dim_t i_Sol=ic+mat->col_block_size*l_col;
1000     for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1001     const dim_t i_Eq=ir+mat->row_block_size*l_row;
1002     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())];
1003     }
1004     }
1005     break;
1006     }
1007     }
1008     } else {
1009     for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1010     if (col_coupleBlock_index[k] == i_col - numMyCols) {
1011     for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1012     const dim_t i_Sol=ic+mat->col_block_size*l_col;
1013     for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1014     const dim_t i_Eq=ir+mat->row_block_size*l_row;
1015     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())];
1016     }
1017     }
1018     break;
1019     }
1020     }
1021 caltinay 3756 }
1022     }
1023     }
1024 caltinay 3760 } else {
1025     for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1026     // across rows of array
1027     for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1028     const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1029     if (i_col < numMyCols) {
1030     for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1031     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1032     {
1033     if (row_coupleBlock_index[k] == i_col) {
1034     for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1035     const dim_t i_Sol=ic+mat->col_block_size*l_col;
1036     for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1037     const dim_t i_Eq=ir+mat->row_block_size*l_row;
1038     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())];
1039     }
1040     }
1041     break;
1042     }
1043     }
1044 caltinay 3756 }
1045     }
1046     }
1047     }
1048     }
1049     }
1050     }
1051    
1052 caltinay 3691 //
1053     // the methods that follow have to be implemented by the subclasses
1054     //
1055    
1056     string RipleyDomain::getDescription() const
1057 caltinay 3670 {
1058 caltinay 3691 throw RipleyException("getDescription() not implemented");
1059 caltinay 3670 }
1060    
1061 caltinay 3691 void RipleyDomain::write(const string& filename) const
1062 caltinay 3670 {
1063 caltinay 3691 throw RipleyException("write() not implemented");
1064 caltinay 3670 }
1065    
1066 caltinay 3691 void RipleyDomain::dump(const string& filename) const
1067 caltinay 3670 {
1068 caltinay 3691 throw RipleyException("dump() not implemented");
1069 caltinay 3670 }
1070    
1071 caltinay 3691 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
1072     {
1073     throw RipleyException("borrowSampleReferenceIDs() not implemented");
1074     }
1075 caltinay 3670
1076 caltinay 3691 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1077 caltinay 3670 {
1078 caltinay 3691 throw RipleyException("interpolateACross() not implemented");
1079 caltinay 3670 }
1080    
1081 caltinay 3691 bool RipleyDomain::probeInterpolationACross(int fsType_source,
1082     const escript::AbstractDomain&, int fsType_target) const
1083 caltinay 3670 {
1084 caltinay 3691 throw RipleyException("probeInterpolationACross() not implemented");
1085 caltinay 3670 }
1086    
1087 caltinay 3691 void RipleyDomain::setToNormal(escript::Data& normal) const
1088 caltinay 3670 {
1089 caltinay 3691 throw RipleyException("setToNormal() not implemented");
1090 caltinay 3670 }
1091    
1092 caltinay 3691 void RipleyDomain::setToSize(escript::Data& size) const
1093 caltinay 3670 {
1094 caltinay 3691 throw RipleyException("setToSize() not implemented");
1095 caltinay 3670 }
1096    
1097 caltinay 3691 bool RipleyDomain::ownSample(int fsType, index_t id) const
1098 caltinay 3670 {
1099 caltinay 3691 throw RipleyException("ownSample() not implemented");
1100 caltinay 3670 }
1101    
1102 caltinay 3691 void RipleyDomain::addPDEToTransportProblem(
1103     escript::AbstractTransportProblem& tp,
1104     escript::Data& source, const escript::Data& M,
1105     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1106     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1107     const escript::Data& d, const escript::Data& y,
1108     const escript::Data& d_contact, const escript::Data& y_contact,
1109     const escript::Data& d_dirac, const escript::Data& y_dirac) const
1110 caltinay 3670 {
1111 caltinay 3691 throw RipleyException("addPDEToTransportProblem() not implemented");
1112 caltinay 3670 }
1113    
1114 caltinay 3691 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
1115     const int blocksize, const escript::FunctionSpace& functionspace,
1116     const int type) const
1117 caltinay 3670 {
1118 caltinay 3691 throw RipleyException("newTransportProblem() not implemented");
1119 caltinay 3670 }
1120    
1121 caltinay 3691 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1122     bool reducedColOrder) const
1123 caltinay 3670 {
1124 caltinay 3691 throw RipleyException("getPattern() not implemented");
1125 caltinay 3670 }
1126    
1127 caltinay 3691 dim_t RipleyDomain::getNumDataPointsGlobal() const
1128 caltinay 3670 {
1129 caltinay 3691 throw RipleyException("getNumDataPointsGlobal() not implemented");
1130 caltinay 3670 }
1131    
1132 caltinay 3697 IndexVector RipleyDomain::getNumNodesPerDim() const
1133     {
1134     throw RipleyException("getNumNodesPerDim() not implemented");
1135     }
1136    
1137     IndexVector RipleyDomain::getNumElementsPerDim() const
1138     {
1139     throw RipleyException("getNumElementsPerDim() not implemented");
1140     }
1141    
1142     IndexVector RipleyDomain::getNumFacesPerBoundary() const
1143     {
1144     throw RipleyException("getNumFacesPerBoundary() not implemented");
1145     }
1146    
1147     IndexVector RipleyDomain::getNodeDistribution() const
1148     {
1149     throw RipleyException("getNodeDistribution() not implemented");
1150     }
1151    
1152 caltinay 3766 IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1153     {
1154     throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1155     }
1156    
1157 caltinay 3697 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1158     {
1159     throw RipleyException("getFirstCoordAndSpacing() not implemented");
1160     }
1161    
1162 caltinay 3691 dim_t RipleyDomain::getNumFaceElements() const
1163 caltinay 3670 {
1164 caltinay 3691 throw RipleyException("getNumFaceElements() not implemented");
1165 caltinay 3670 }
1166    
1167 caltinay 3691 dim_t RipleyDomain::getNumElements() const
1168 caltinay 3670 {
1169 caltinay 3691 throw RipleyException("getNumElements() not implemented");
1170 caltinay 3670 }
1171    
1172 caltinay 3691 dim_t RipleyDomain::getNumNodes() const
1173 caltinay 3670 {
1174 caltinay 3691 throw RipleyException("getNumNodes() not implemented");
1175 caltinay 3670 }
1176    
1177 caltinay 3750 dim_t RipleyDomain::getNumDOF() const
1178     {
1179     throw RipleyException("getNumDOF() not implemented");
1180     }
1181    
1182 caltinay 3756 dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1183     {
1184     throw RipleyException("insertNeighbourNodes() not implemented");
1185     }
1186    
1187 caltinay 3691 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1188 caltinay 3670 {
1189 caltinay 3691 throw RipleyException("assembleCoordinates() not implemented");
1190 caltinay 3670 }
1191    
1192 caltinay 3764 void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1193     {
1194     throw RipleyException("assembleGradient() not implemented");
1195     }
1196    
1197     void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1198     {
1199     throw RipleyException("assembleIntegrate() not implemented");
1200     }
1201    
1202 caltinay 3748 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1203     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1204 caltinay 3769 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1205 caltinay 3748 {
1206     throw RipleyException("assemblePDESingle() not implemented");
1207     }
1208    
1209 caltinay 3769 void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1210 caltinay 3776 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1211 caltinay 3769 {
1212     throw RipleyException("assemblePDEBoundarySingle() not implemented");
1213     }
1214    
1215 caltinay 3761 void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1216     escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1217     const escript::Data& C, const escript::Data& D, const escript::Data& X,
1218 caltinay 3769 const escript::Data& Y) const
1219 caltinay 3761 {
1220     throw RipleyException("assemblePDESingleReduced() not implemented");
1221     }
1222    
1223 caltinay 3769 void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1224 caltinay 3776 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1225 caltinay 3769 {
1226     throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1227     }
1228    
1229 caltinay 3748 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1230     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1231 caltinay 3769 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1232 caltinay 3748 {
1233     throw RipleyException("assemblePDESystem() not implemented");
1234     }
1235    
1236 caltinay 3769 void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1237 caltinay 3776 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1238 caltinay 3769 {
1239     throw RipleyException("assemblePDEBoundarySystem() not implemented");
1240     }
1241    
1242 caltinay 3761 void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1243     escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1244     const escript::Data& C, const escript::Data& D, const escript::Data& X,
1245 caltinay 3769 const escript::Data& Y) const
1246 caltinay 3761 {
1247     throw RipleyException("assemblePDESystemReduced() not implemented");
1248     }
1249    
1250 caltinay 3769 void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1251 caltinay 3776 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1252 caltinay 3769 {
1253     throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1254     }
1255    
1256 caltinay 3711 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1257 caltinay 3701 {
1258     throw RipleyException("interpolateNodesOnElements() not implemented");
1259     }
1260 caltinay 3691
1261 caltinay 3711 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1262 caltinay 3701 {
1263     throw RipleyException("interpolateNodesOnFaces() not implemented");
1264     }
1265    
1266 caltinay 3750 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1267     {
1268     throw RipleyException("nodesToDOF() not implemented");
1269     }
1270 caltinay 3701
1271 caltinay 3754 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1272     {
1273     throw RipleyException("dofToNodes() not implemented");
1274     }
1275    
1276 caltinay 3670 } // end of namespace ripley
1277    

  ViewVC Help
Powered by ViewVC 1.1.26