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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3791 - (hide annotations)
Wed Feb 1 05:10:22 2012 UTC (7 years, 7 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 45625 byte(s)
Getting ready to merge with trunk.

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 3791 msg << "getDataShape: Invalid function space type " << fsType
124 caltinay 3740 << " 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 caltinay 3791 msg << "probeInterpolationOnDomain: Invalid function space 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 caltinay 3791 msg << "probeInterpolationOnDomain: Invalid function space 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 3791 throw RipleyException("interpolateOnDomain: Cannot interpolate reduced data to non-reduced data.");
298 caltinay 3744 } 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 caltinay 3791 msg << "setToGradient: not supported for "
445 caltinay 3764 << 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 caltinay 3791 throw RipleyException("setToIntegrals: illegal domain of integration kernel");
471 caltinay 3764
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 caltinay 3791 msg << "setToIntegrals: not supported for "
491 caltinay 3764 << 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 3791 msg << "isCellOriented: invalid function space type " << fsType
517 caltinay 3740 << " 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 3791 msg << "canTag: invalid function space type " << fsType << " on "
540 caltinay 3740 << 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 3791 msg << "setTags: invalid function space type " << fsType;
567 caltinay 3722 throw RipleyException(msg.str());
568     }
569     }
570     if (target->size() != num) {
571     target->assign(num, -1);
572     }
573    
574     escript::Data& mask=*const_cast<escript::Data*>(&cMask);
575     #pragma omp parallel for
576     for (index_t i=0; i<num; i++) {
577     if (mask.getSampleDataRO(i)[0] > 0) {
578     (*target)[i]=newTag;
579     }
580     }
581     updateTagsInUse(fsType);
582     }
583    
584     int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
585     {
586     switch(fsType) {
587     case Nodes:
588     if (m_nodeTags.size() > sampleNo)
589     return m_nodeTags[sampleNo];
590     break;
591     case Elements:
592     case ReducedElements:
593     if (m_elementTags.size() > sampleNo)
594     return m_elementTags[sampleNo];
595     break;
596     case FaceElements:
597     case ReducedFaceElements:
598     if (m_faceTags.size() > sampleNo)
599     return m_faceTags[sampleNo];
600     break;
601     default: {
602     stringstream msg;
603 caltinay 3791 msg << "getTagFromSampleNo: invalid function space type " << fsType;
604 caltinay 3722 throw RipleyException(msg.str());
605     }
606     }
607     return -1;
608     }
609    
610     int RipleyDomain::getNumberOfTagsInUse(int fsType) const
611     {
612     switch(fsType) {
613     case Nodes:
614     return m_nodeTagsInUse.size();
615     case Elements:
616     case ReducedElements:
617     return m_elementTagsInUse.size();
618     case FaceElements:
619     case ReducedFaceElements:
620     return m_faceTagsInUse.size();
621     default: {
622     stringstream msg;
623 caltinay 3791 msg << "getNumberOfTagsInUse: invalid function space type "
624     << fsType;
625 caltinay 3722 throw RipleyException(msg.str());
626     }
627     }
628     }
629    
630     const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
631     {
632     switch(fsType) {
633     case Nodes:
634     return &m_nodeTagsInUse[0];
635     case Elements:
636     case ReducedElements:
637     return &m_elementTagsInUse[0];
638     case FaceElements:
639     case ReducedFaceElements:
640     return &m_faceTagsInUse[0];
641     default: {
642     stringstream msg;
643 caltinay 3791 msg << "borrowListOfTagsInUse: invalid function space type "
644     << fsType;
645 caltinay 3722 throw RipleyException(msg.str());
646     }
647     }
648     }
649    
650 caltinay 3691 void RipleyDomain::Print_Mesh_Info(const bool full) const
651 caltinay 3670 {
652 caltinay 3691 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
653     << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
654     cout << "Number of dimensions: " << m_numDim << endl;
655 caltinay 3670
656 caltinay 3691 // write tags
657     if (m_tagMap.size() > 0) {
658     cout << "Tags:" << endl;
659     TagMap::const_iterator it;
660     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
661     cout << " " << setw(5) << it->second << " "
662     << it->first << endl;
663     }
664 caltinay 3670 }
665     }
666    
667     int RipleyDomain::getSystemMatrixTypeId(const int solver,
668     const int preconditioner, const int package, const bool symmetry) const
669     {
670 caltinay 3691 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
671     package, symmetry, m_mpiInfo);
672 caltinay 3670 }
673    
674     int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
675     const int package, const bool symmetry) const
676     {
677 caltinay 3691 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
678     package, symmetry, m_mpiInfo);
679 caltinay 3670 }
680    
681 caltinay 3691 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
682     const escript::FunctionSpace& row_functionspace,
683     const int column_blocksize,
684     const escript::FunctionSpace& column_functionspace,
685     const int type) const
686 caltinay 3670 {
687 caltinay 3691 bool reduceRowOrder=false;
688     bool reduceColOrder=false;
689     // is the domain right?
690     const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
691 caltinay 3764 if (row_domain != *this)
692 caltinay 3791 throw RipleyException("newSystemMatrix: domain of row function space does not match the domain of matrix generator");
693 caltinay 3691 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
694 caltinay 3764 if (col_domain != *this)
695 caltinay 3791 throw RipleyException("newSystemMatrix: domain of column function space does not match the domain of matrix generator");
696 caltinay 3691 // is the function space type right?
697 caltinay 3750 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
698 caltinay 3691 reduceRowOrder=true;
699 caltinay 3750 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
700 caltinay 3791 throw RipleyException("newSystemMatrix: illegal function space type for system matrix rows");
701 caltinay 3750 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
702 caltinay 3691 reduceColOrder=true;
703 caltinay 3750 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
704 caltinay 3791 throw RipleyException("newSystemMatrix: illegal function space type for system matrix columns");
705 caltinay 3691
706     // generate matrix
707     Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
708     Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
709     row_blocksize, column_blocksize, FALSE);
710     paso::checkPasoError();
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 caltinay 3791 throw RipleyException("addPDEToSystem: Ripley does not support contact elements");
726 caltinay 3748
727     paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
728     if (!sma)
729 caltinay 3791 throw RipleyException("addPDEToSystem: Ripley only accepts Paso system matrices");
730 caltinay 3748
731     if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
732 caltinay 3791 throw RipleyException("addPDEToSystem: right hand side coefficients are provided but no right hand side vector given");
733 caltinay 3748
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 caltinay 3791 throw RipleyException("addPDEToSystem: matrix row block size and number of components of right hand side don't match");
756 caltinay 3748
757     const int numEq=S->logical_row_block_size;
758     const int numComp=S->logical_col_block_size;
759    
760     if (numEq != numComp)
761 caltinay 3791 throw RipleyException("addPDEToSystem: number of equations and number of solutions don't match");
762 caltinay 3748 //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 caltinay 3791 throw RipleyException("addPDEToRHS: Ripley does not support contact elements");
789 caltinay 3760
790     if (rhs.isEmpty()) {
791     if (!X.isEmpty() || !Y.isEmpty())
792 caltinay 3791 throw RipleyException("addPDEToRHS: right hand side coefficients are provided but no right hand side vector given");
793 caltinay 3760 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 caltinay 3791 throw RipleyException("setNewX(): operation not supported");
809 caltinay 3697 }
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::addPDEToTransportProblem(
1101     escript::AbstractTransportProblem& tp,
1102     escript::Data& source, const escript::Data& M,
1103     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1104     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1105     const escript::Data& d, const escript::Data& y,
1106     const escript::Data& d_contact, const escript::Data& y_contact,
1107     const escript::Data& d_dirac, const escript::Data& y_dirac) const
1108 caltinay 3670 {
1109 caltinay 3691 throw RipleyException("addPDEToTransportProblem() not implemented");
1110 caltinay 3670 }
1111    
1112 caltinay 3691 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
1113     const int blocksize, const escript::FunctionSpace& functionspace,
1114     const int type) const
1115 caltinay 3670 {
1116 caltinay 3691 throw RipleyException("newTransportProblem() not implemented");
1117 caltinay 3670 }
1118    
1119 caltinay 3691 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1120     bool reducedColOrder) const
1121 caltinay 3670 {
1122 caltinay 3691 throw RipleyException("getPattern() not implemented");
1123 caltinay 3670 }
1124    
1125 caltinay 3691 dim_t RipleyDomain::getNumDataPointsGlobal() const
1126 caltinay 3670 {
1127 caltinay 3691 throw RipleyException("getNumDataPointsGlobal() not implemented");
1128 caltinay 3670 }
1129    
1130 caltinay 3697 IndexVector RipleyDomain::getNumNodesPerDim() const
1131     {
1132     throw RipleyException("getNumNodesPerDim() not implemented");
1133     }
1134    
1135     IndexVector RipleyDomain::getNumElementsPerDim() const
1136     {
1137     throw RipleyException("getNumElementsPerDim() not implemented");
1138     }
1139    
1140     IndexVector RipleyDomain::getNumFacesPerBoundary() const
1141     {
1142     throw RipleyException("getNumFacesPerBoundary() not implemented");
1143     }
1144    
1145     IndexVector RipleyDomain::getNodeDistribution() const
1146     {
1147     throw RipleyException("getNodeDistribution() not implemented");
1148     }
1149    
1150 caltinay 3766 IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1151     {
1152     throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1153     }
1154    
1155 caltinay 3697 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1156     {
1157     throw RipleyException("getFirstCoordAndSpacing() not implemented");
1158     }
1159    
1160 caltinay 3691 dim_t RipleyDomain::getNumFaceElements() const
1161 caltinay 3670 {
1162 caltinay 3691 throw RipleyException("getNumFaceElements() not implemented");
1163 caltinay 3670 }
1164    
1165 caltinay 3691 dim_t RipleyDomain::getNumElements() const
1166 caltinay 3670 {
1167 caltinay 3691 throw RipleyException("getNumElements() not implemented");
1168 caltinay 3670 }
1169    
1170 caltinay 3691 dim_t RipleyDomain::getNumNodes() const
1171 caltinay 3670 {
1172 caltinay 3691 throw RipleyException("getNumNodes() not implemented");
1173 caltinay 3670 }
1174    
1175 caltinay 3750 dim_t RipleyDomain::getNumDOF() const
1176     {
1177     throw RipleyException("getNumDOF() not implemented");
1178     }
1179    
1180 caltinay 3756 dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1181     {
1182     throw RipleyException("insertNeighbourNodes() not implemented");
1183     }
1184    
1185 caltinay 3691 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1186 caltinay 3670 {
1187 caltinay 3691 throw RipleyException("assembleCoordinates() not implemented");
1188 caltinay 3670 }
1189    
1190 caltinay 3764 void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1191     {
1192     throw RipleyException("assembleGradient() not implemented");
1193     }
1194    
1195     void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1196     {
1197     throw RipleyException("assembleIntegrate() not implemented");
1198     }
1199    
1200 caltinay 3748 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1201     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1202 caltinay 3769 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1203 caltinay 3748 {
1204     throw RipleyException("assemblePDESingle() not implemented");
1205     }
1206    
1207 caltinay 3769 void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1208 caltinay 3776 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1209 caltinay 3769 {
1210     throw RipleyException("assemblePDEBoundarySingle() not implemented");
1211     }
1212    
1213 caltinay 3761 void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1214     escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1215     const escript::Data& C, const escript::Data& D, const escript::Data& X,
1216 caltinay 3769 const escript::Data& Y) const
1217 caltinay 3761 {
1218     throw RipleyException("assemblePDESingleReduced() not implemented");
1219     }
1220    
1221 caltinay 3769 void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1222 caltinay 3776 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1223 caltinay 3769 {
1224     throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1225     }
1226    
1227 caltinay 3748 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1228     const escript::Data& A, const escript::Data& B, const escript::Data& C,
1229 caltinay 3769 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1230 caltinay 3748 {
1231     throw RipleyException("assemblePDESystem() not implemented");
1232     }
1233    
1234 caltinay 3769 void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1235 caltinay 3776 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1236 caltinay 3769 {
1237     throw RipleyException("assemblePDEBoundarySystem() not implemented");
1238     }
1239    
1240 caltinay 3761 void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1241     escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1242     const escript::Data& C, const escript::Data& D, const escript::Data& X,
1243 caltinay 3769 const escript::Data& Y) const
1244 caltinay 3761 {
1245     throw RipleyException("assemblePDESystemReduced() not implemented");
1246     }
1247    
1248 caltinay 3769 void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1249 caltinay 3776 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1250 caltinay 3769 {
1251     throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1252     }
1253    
1254 caltinay 3711 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1255 caltinay 3701 {
1256     throw RipleyException("interpolateNodesOnElements() not implemented");
1257     }
1258 caltinay 3691
1259 caltinay 3711 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1260 caltinay 3701 {
1261     throw RipleyException("interpolateNodesOnFaces() not implemented");
1262     }
1263    
1264 caltinay 3750 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1265     {
1266     throw RipleyException("nodesToDOF() not implemented");
1267     }
1268 caltinay 3701
1269 caltinay 3754 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1270     {
1271     throw RipleyException("dofToNodes() not implemented");
1272     }
1273    
1274 caltinay 3670 } // end of namespace ripley
1275    

  ViewVC Help
Powered by ViewVC 1.1.26