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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3761 - (hide annotations)
Mon Jan 9 06:50:33 2012 UTC (7 years, 8 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 42425 byte(s)
Can solve full order systems of PDEs in 2D and 3D now.

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

  ViewVC Help
Powered by ViewVC 1.1.26