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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3780 - (hide annotations)
Mon Jan 23 06:26:40 2012 UTC (7 years, 8 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 45825 byte(s)
Removed addPDEToLumpedSystem() as it is handled at python level.
Ensured tests for setX() don't fail if unsupported.

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

  ViewVC Help
Powered by ViewVC 1.1.26