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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3748 - (hide annotations)
Thu Dec 15 07:36:19 2011 UTC (7 years, 9 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 28471 byte(s)
PDE assembly in serial and 2D seems to be doing what it's supposed to when
boundary elements are not involved.

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 3670 case Nodes:
65     case ReducedNodes:
66     case Elements:
67     case ReducedElements:
68     case FaceElements:
69     case ReducedFaceElements:
70     case Points:
71 caltinay 3691 return true;
72     default:
73 caltinay 3670 break;
74     }
75 caltinay 3691 return false;
76 caltinay 3670 }
77    
78 caltinay 3691 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
79 caltinay 3670 {
80 caltinay 3691 switch (fsType) {
81     case Nodes: return "Ripley_Nodes";
82     case ReducedNodes: return "Ripley_Reduced_Nodes";
83     case Elements: return "Ripley_Elements";
84     case ReducedElements: return "Ripley_Reduced_Elements";
85     case FaceElements: return "Ripley_Face_Elements";
86     case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";
87     case Points: return "Ripley_Points";
88     default:
89 caltinay 3670 break;
90     }
91 caltinay 3691 return "Invalid function space type code";
92 caltinay 3670 }
93    
94 caltinay 3691 pair<int,int> RipleyDomain::getDataShape(int fsType) const
95 caltinay 3670 {
96 caltinay 3691 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
97     switch (fsType) {
98 caltinay 3670 case Nodes:
99 caltinay 3744 case ReducedNodes: //FIXME: reduced
100     return pair<int,int>(1, getNumNodes());
101 caltinay 3670 case Elements:
102 caltinay 3691 return pair<int,int>(ptsPerSample, getNumElements());
103 caltinay 3670 case FaceElements:
104 caltinay 3691 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
105 caltinay 3711 case ReducedElements:
106     return pair<int,int>(1, getNumElements());
107     case ReducedFaceElements:
108     return pair<int,int>(1, getNumFaceElements());
109 caltinay 3697 case Points:
110 caltinay 3733 return pair<int,int>(1, 0); //FIXME: dirac
111 caltinay 3670 default:
112     break;
113     }
114    
115 caltinay 3691 stringstream msg;
116 caltinay 3740 msg << "getDataShape(): Unsupported function space type " << fsType
117     << " for " << getDescription();
118 caltinay 3691 throw RipleyException(msg.str());
119 caltinay 3670 }
120    
121 caltinay 3691 string RipleyDomain::showTagNames() const
122 caltinay 3670 {
123 caltinay 3691 stringstream ret;
124     TagMap::const_iterator it;
125     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
126     if (it!=m_tagMap.begin()) ret << ", ";
127     ret << it->first;
128 caltinay 3670 }
129 caltinay 3691 return ret.str();
130 caltinay 3670 }
131    
132     bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
133     {
134 caltinay 3691 /*
135     The idea is to use equivalence classes (i.e. types which can be
136     interpolated back and forth):
137 caltinay 3746 class 0: Nodes
138     class 1: ReducedNodes
139 caltinay 3740 class 2: Points
140     class 3: Elements
141     class 4: ReducedElements
142     class 5: FaceElements
143     class 6: ReducedFaceElements
144 caltinay 3670
145 caltinay 3691 There is also a set of lines. Interpolation is possible down a line but not
146     between lines.
147 caltinay 3740 class 0 and 1 belong to all lines so aren't considered.
148     line 0: class 2
149     line 1: class 3,4
150     line 2: class 5,6
151 caltinay 3670 */
152     if (fs.empty())
153     return false;
154 caltinay 3740 vector<bool> hasclass(7, false);
155     vector<int> hasline(3, 0);
156     for (size_t i=0; i<fs.size(); ++i) {
157 caltinay 3670 switch (fs[i]) {
158 caltinay 3746 case Nodes:
159 caltinay 3740 hasclass[0]=true;
160 caltinay 3670 break;
161 caltinay 3746 case ReducedNodes:
162 caltinay 3740 hasclass[1]=true;
163 caltinay 3670 break;
164     case Points:
165     hasline[0]=1;
166 caltinay 3740 hasclass[2]=true;
167 caltinay 3670 break;
168     case Elements:
169     hasline[1]=1;
170     break;
171     case ReducedElements:
172 caltinay 3740 hasclass[4]=true;
173 caltinay 3670 hasline[1]=1;
174     break;
175     case FaceElements:
176     hasline[2]=1;
177     break;
178     case ReducedFaceElements:
179 caltinay 3740 hasclass[6]=true;
180 caltinay 3670 hasline[2]=1;
181     break;
182     default:
183     return false;
184     }
185     }
186 caltinay 3740 int numLines=hasline[0]+hasline[1]+hasline[2];
187 caltinay 3670
188     // fail if we have more than one leaf group
189     // = there are at least two branches we can't interpolate between
190 caltinay 3740 if (numLines > 1)
191 caltinay 3670 return false;
192 caltinay 3740 else if (numLines==1) {
193 caltinay 3691 if (hasline[0]==1)
194 caltinay 3670 resultcode=Points;
195 caltinay 3691 else if (hasline[1]==1) {
196 caltinay 3740 if (hasclass[4])
197 caltinay 3670 resultcode=ReducedElements;
198 caltinay 3691 else
199 caltinay 3670 resultcode=Elements;
200 caltinay 3740 } else { // hasline[2]==1
201     if (hasclass[6])
202 caltinay 3670 resultcode=ReducedFaceElements;
203 caltinay 3691 else
204 caltinay 3670 resultcode=FaceElements;
205 caltinay 3740 }
206     } else { // numLines==0
207     if (hasclass[1])
208 caltinay 3746 resultcode=ReducedNodes;
209     else
210     resultcode=Nodes;
211 caltinay 3670 }
212     return true;
213     }
214    
215 caltinay 3732 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
216     int fsType_target) const
217     {
218     if (fsType_target != Nodes &&
219     fsType_target != ReducedNodes &&
220     fsType_target != Elements &&
221     fsType_target != ReducedElements &&
222     fsType_target != FaceElements &&
223     fsType_target != ReducedFaceElements &&
224     fsType_target != Points) {
225     stringstream msg;
226     msg << "probeInterpolationOnDomain(): Invalid functionspace type "
227 caltinay 3740 << fsType_target << " for " << getDescription();
228 caltinay 3732 throw RipleyException(msg.str());
229     }
230    
231     switch (fsType_source) {
232     case Nodes:
233     return true;
234     case ReducedNodes:
235 caltinay 3746 return (fsType_target != Nodes);
236 caltinay 3732 case Elements:
237     return (fsType_target==Elements ||
238     fsType_target==ReducedElements);
239     case ReducedElements:
240     return (fsType_target==ReducedElements);
241     case FaceElements:
242     return (fsType_target==FaceElements ||
243     fsType_target==ReducedFaceElements);
244     case ReducedFaceElements:
245     return (fsType_target==ReducedFaceElements);
246     case Points:
247     return (fsType_target==Points);
248    
249     default: {
250     stringstream msg;
251     msg << "probeInterpolationOnDomain(): Invalid functionspace type "
252 caltinay 3740 << fsType_source << " for " << getDescription();
253 caltinay 3732 throw RipleyException(msg.str());
254     }
255     }
256     }
257    
258 caltinay 3701 void RipleyDomain::interpolateOnDomain(escript::Data& target,
259     const escript::Data& in) const
260     {
261     const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
262     const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
263     if (inDomain != *this)
264     throw RipleyException("Illegal domain of interpolant");
265     if (targetDomain != *this)
266     throw RipleyException("Illegal domain of interpolation target");
267    
268     stringstream msg;
269     msg << "interpolateOnDomain() not implemented for function space "
270     << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
271     << " -> "
272     << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
273    
274 caltinay 3744 const int inFS = in.getFunctionSpace().getTypeCode();
275     const int outFS = target.getFunctionSpace().getTypeCode();
276 caltinay 3701
277 caltinay 3744 // simplest case: 1:1 copy
278     if (inFS==outFS) {
279     copyData(target, *const_cast<escript::Data*>(&in));
280     // not allowed: reduced->non-reduced
281 caltinay 3746 } else if (inFS==ReducedNodes && outFS==Nodes) {
282 caltinay 3744 throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
283     } else if ((inFS==Elements && outFS==ReducedElements)
284     || (inFS==FaceElements && outFS==ReducedFaceElements)) {
285     averageData(target, *const_cast<escript::Data*>(&in));
286     } else {
287     switch (inFS) {
288     case Nodes:
289 caltinay 3746 case ReducedNodes:
290 caltinay 3744 switch (outFS) {
291     case Nodes:
292     case ReducedNodes: //FIXME: reduced
293 caltinay 3746 copyData(target, *const_cast<escript::Data*>(&in));
294 caltinay 3744 break;
295 caltinay 3701
296 caltinay 3744 case Elements:
297     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
298     break;
299 caltinay 3711
300 caltinay 3744 case ReducedElements:
301     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
302     break;
303 caltinay 3701
304 caltinay 3744 case FaceElements:
305     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
306     break;
307 caltinay 3711
308 caltinay 3744 case ReducedFaceElements:
309     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
310     break;
311    
312     default:
313     throw RipleyException(msg.str());
314     }
315     break;
316     default:
317     throw RipleyException(msg.str());
318     }
319 caltinay 3701 }
320     }
321    
322 caltinay 3691 escript::Data RipleyDomain::getX() const
323 caltinay 3670 {
324 caltinay 3691 return escript::continuousFunction(*this).getX();
325     }
326    
327     escript::Data RipleyDomain::getNormal() const
328     {
329     return escript::functionOnBoundary(*this).getNormal();
330     }
331    
332     escript::Data RipleyDomain::getSize() const
333     {
334     return escript::function(*this).getSize();
335     }
336    
337     void RipleyDomain::setToX(escript::Data& arg) const
338     {
339     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
340     *(arg.getFunctionSpace().getDomain()));
341     if (argDomain != *this)
342     throw RipleyException("setToX: Illegal domain of data point locations");
343     if (!arg.isExpanded())
344     throw RipleyException("setToX: Expanded Data object expected");
345    
346     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
347     assembleCoordinates(arg);
348     } else {
349     // interpolate the result
350     escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
351     assembleCoordinates(contData);
352     interpolateOnDomain(arg, contData);
353 caltinay 3670 }
354 caltinay 3691 }
355 caltinay 3670
356 caltinay 3691 bool RipleyDomain::isCellOriented(int fsType) const
357     {
358     switch(fsType) {
359 caltinay 3670 case Nodes:
360 caltinay 3691 return false;
361 caltinay 3670 case Elements:
362 caltinay 3691 case FaceElements:
363     case Points:
364 caltinay 3670 case ReducedElements:
365     case ReducedFaceElements:
366 caltinay 3691 return true;
367 caltinay 3670 default:
368 caltinay 3691 break;
369 caltinay 3670 }
370 caltinay 3691 stringstream msg;
371 caltinay 3740 msg << "isCellOriented(): Illegal function space type " << fsType
372     << " on " << getDescription();
373 caltinay 3691 throw RipleyException(msg.str());
374 caltinay 3670 }
375    
376 caltinay 3722 bool RipleyDomain::canTag(int fsType) const
377     {
378     switch(fsType) {
379     case Nodes:
380     case Elements:
381     case ReducedElements:
382     case FaceElements:
383     case ReducedFaceElements:
384     return true;
385     case Points:
386     return false;
387     default:
388     break;
389     }
390     stringstream msg;
391 caltinay 3740 msg << "canTag(): Illegal function space type " << fsType << " on "
392     << getDescription();
393 caltinay 3722 throw RipleyException(msg.str());
394     }
395    
396     void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
397     {
398     IndexVector* target=NULL;
399     dim_t num=0;
400    
401     switch(fsType) {
402     case Nodes:
403     num=getNumNodes();
404     target=&m_nodeTags;
405     break;
406     case Elements:
407     case ReducedElements:
408     num=getNumElements();
409     target=&m_elementTags;
410     break;
411     case FaceElements:
412     case ReducedFaceElements:
413     num=getNumFaceElements();
414     target=&m_faceTags;
415     break;
416     default: {
417     stringstream msg;
418 caltinay 3740 msg << "setTags(): not implemented for "
419     << functionSpaceTypeAsString(fsType);
420 caltinay 3722 throw RipleyException(msg.str());
421     }
422     }
423     if (target->size() != num) {
424     target->assign(num, -1);
425     }
426    
427     escript::Data& mask=*const_cast<escript::Data*>(&cMask);
428     #pragma omp parallel for
429     for (index_t i=0; i<num; i++) {
430     if (mask.getSampleDataRO(i)[0] > 0) {
431     (*target)[i]=newTag;
432     }
433     }
434     updateTagsInUse(fsType);
435     }
436    
437     int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
438     {
439     switch(fsType) {
440     case Nodes:
441     if (m_nodeTags.size() > sampleNo)
442     return m_nodeTags[sampleNo];
443     break;
444     case Elements:
445     case ReducedElements:
446     if (m_elementTags.size() > sampleNo)
447     return m_elementTags[sampleNo];
448     break;
449     case FaceElements:
450     case ReducedFaceElements:
451     if (m_faceTags.size() > sampleNo)
452     return m_faceTags[sampleNo];
453     break;
454     default: {
455     stringstream msg;
456 caltinay 3740 msg << "getTagFromSampleNo(): not implemented for "
457     << functionSpaceTypeAsString(fsType);
458 caltinay 3722 throw RipleyException(msg.str());
459     }
460     }
461     return -1;
462     }
463    
464     int RipleyDomain::getNumberOfTagsInUse(int fsType) const
465     {
466     switch(fsType) {
467     case Nodes:
468     return m_nodeTagsInUse.size();
469     case Elements:
470     case ReducedElements:
471     return m_elementTagsInUse.size();
472     case FaceElements:
473     case ReducedFaceElements:
474     return m_faceTagsInUse.size();
475     default: {
476     stringstream msg;
477 caltinay 3740 msg << "getNumberOfTagsInUse(): not implemented for "
478     << functionSpaceTypeAsString(fsType);
479 caltinay 3722 throw RipleyException(msg.str());
480     }
481     }
482     }
483    
484     const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
485     {
486     switch(fsType) {
487     case Nodes:
488     return &m_nodeTagsInUse[0];
489     case Elements:
490     case ReducedElements:
491     return &m_elementTagsInUse[0];
492     case FaceElements:
493     case ReducedFaceElements:
494     return &m_faceTagsInUse[0];
495     default: {
496     stringstream msg;
497 caltinay 3740 msg << "borrowListOfTagsInUse(): not implemented for "
498     << functionSpaceTypeAsString(fsType);
499 caltinay 3722 throw RipleyException(msg.str());
500     }
501     }
502     }
503    
504 caltinay 3691 void RipleyDomain::Print_Mesh_Info(const bool full) const
505 caltinay 3670 {
506 caltinay 3691 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
507     << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
508     cout << "Number of dimensions: " << m_numDim << endl;
509 caltinay 3670
510 caltinay 3691 // write tags
511     if (m_tagMap.size() > 0) {
512     cout << "Tags:" << endl;
513     TagMap::const_iterator it;
514     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
515     cout << " " << setw(5) << it->second << " "
516     << it->first << endl;
517     }
518 caltinay 3670 }
519     }
520    
521     int RipleyDomain::getSystemMatrixTypeId(const int solver,
522     const int preconditioner, const int package, const bool symmetry) const
523     {
524 caltinay 3691 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
525     package, symmetry, m_mpiInfo);
526 caltinay 3670 }
527    
528     int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
529     const int package, const bool symmetry) const
530     {
531 caltinay 3691 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
532     package, symmetry, m_mpiInfo);
533 caltinay 3670 }
534    
535 caltinay 3691 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
536     const escript::FunctionSpace& row_functionspace,
537     const int column_blocksize,
538     const escript::FunctionSpace& column_functionspace,
539     const int type) const
540 caltinay 3670 {
541 caltinay 3691 bool reduceRowOrder=false;
542     bool reduceColOrder=false;
543     // is the domain right?
544     const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
545     if (row_domain!=*this)
546 caltinay 3740 throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
547 caltinay 3691 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
548     if (col_domain!=*this)
549 caltinay 3740 throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
550 caltinay 3691 // is the function space type right?
551 caltinay 3746 if (row_functionspace.getTypeCode()==ReducedNodes)
552 caltinay 3691 reduceRowOrder=true;
553 caltinay 3746 else if (row_functionspace.getTypeCode()!=Nodes)
554 caltinay 3740 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
555 caltinay 3746 if (column_functionspace.getTypeCode()==ReducedNodes)
556 caltinay 3691 reduceColOrder=true;
557 caltinay 3746 else if (column_functionspace.getTypeCode()!=Nodes)
558 caltinay 3740 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
559 caltinay 3691
560     // generate matrix
561     Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
562     Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
563     row_blocksize, column_blocksize, FALSE);
564     paso::checkPasoError();
565     Paso_SystemMatrixPattern_free(pattern);
566     escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
567     row_functionspace, column_blocksize, column_functionspace));
568     return sma;
569 caltinay 3670 }
570    
571 caltinay 3748 void RipleyDomain::addPDEToSystem(
572     escript::AbstractSystemMatrix& mat, escript::Data& rhs,
573     const escript::Data& A, const escript::Data& B, const escript::Data& C,
574     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
575     const escript::Data& d, const escript::Data& y,
576     const escript::Data& d_contact, const escript::Data& y_contact,
577     const escript::Data& d_dirac,const escript::Data& y_dirac) const
578     {
579     if (!d_contact.isEmpty() || !y_contact.isEmpty())
580     throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
581    
582     paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
583     if (!sma)
584     throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
585    
586     if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
587     throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
588    
589     //TODO: more input param checks (shape, function space etc)
590    
591     Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
592    
593     if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
594     throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
595    
596     const int numEq=S->logical_row_block_size;
597     const int numComp=S->logical_col_block_size;
598    
599     if (numEq != numComp)
600     throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
601     //TODO: more system matrix checks
602    
603     if (numEq==1)
604     assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
605     else
606     assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
607     }
608    
609 caltinay 3697 void RipleyDomain::setNewX(const escript::Data& arg)
610     {
611     throw RipleyException("setNewX(): Operation not supported");
612     }
613    
614 caltinay 3722 //protected
615 caltinay 3744 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
616 caltinay 3722 {
617     const dim_t numComp = in.getDataPointSize();
618 caltinay 3744 const dim_t dpp = in.getNumDataPointsPerSample();
619 caltinay 3740 out.requireWrite();
620 caltinay 3722 #pragma omp parallel for
621     for (index_t i=0; i<in.getNumSamples(); i++) {
622     const double* src = in.getSampleDataRO(i);
623 caltinay 3744 double* dest = out.getSampleDataRW(i);
624     for (index_t c=0; c<numComp; c++) {
625     double res=0.;
626     for (index_t q=0; q<dpp; q++)
627     res+=src[c+q*numComp];
628     *dest++ = res/dpp;
629     }
630     }
631     }
632    
633     //protected
634     void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
635     {
636     const dim_t numComp = in.getDataPointSize();
637     out.requireWrite();
638     #pragma omp parallel for
639     for (index_t i=0; i<in.getNumSamples(); i++) {
640     const double* src = in.getSampleDataRO(i);
641 caltinay 3722 copy(src, src+numComp, out.getSampleDataRW(i));
642     }
643     }
644    
645     //protected
646     void RipleyDomain::updateTagsInUse(int fsType) const
647     {
648     IndexVector* tagsInUse=NULL;
649     const IndexVector* tags=NULL;
650     switch(fsType) {
651     case Nodes:
652     tags=&m_nodeTags;
653     tagsInUse=&m_nodeTagsInUse;
654     break;
655     case Elements:
656     case ReducedElements:
657     tags=&m_elementTags;
658     tagsInUse=&m_elementTagsInUse;
659     break;
660     case FaceElements:
661     case ReducedFaceElements:
662     tags=&m_faceTags;
663     tagsInUse=&m_faceTagsInUse;
664     break;
665     default:
666     return;
667     }
668    
669     // gather global unique tag values from tags into tagsInUse
670     tagsInUse->clear();
671     index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
672    
673     while (true) {
674     // find smallest value bigger than lastFoundValue
675     minFoundValue = INDEX_T_MAX;
676     #pragma omp parallel private(local_minFoundValue)
677     {
678     local_minFoundValue = minFoundValue;
679     #pragma omp for schedule(static) nowait
680     for (size_t i = 0; i < tags->size(); i++) {
681     const index_t v = (*tags)[i];
682     if ((v > lastFoundValue) && (v < local_minFoundValue))
683     local_minFoundValue = v;
684     }
685     #pragma omp critical
686     {
687     if (local_minFoundValue < minFoundValue)
688     minFoundValue = local_minFoundValue;
689     }
690     }
691     #ifdef ESYS_MPI
692     local_minFoundValue = minFoundValue;
693     MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
694     #endif
695 caltinay 3740
696 caltinay 3722 // if we found a new value add it to the tagsInUse vector
697     if (minFoundValue < INDEX_T_MAX) {
698     tagsInUse->push_back(minFoundValue);
699     lastFoundValue = minFoundValue;
700     } else
701     break;
702     }
703     }
704    
705 caltinay 3691 //
706     // the methods that follow have to be implemented by the subclasses
707     //
708    
709     string RipleyDomain::getDescription() const
710 caltinay 3670 {
711 caltinay 3691 throw RipleyException("getDescription() not implemented");
712 caltinay 3670 }
713    
714 caltinay 3691 void RipleyDomain::write(const string& filename) const
715 caltinay 3670 {
716 caltinay 3691 throw RipleyException("write() not implemented");
717 caltinay 3670 }
718    
719 caltinay 3691 void RipleyDomain::dump(const string& filename) const
720 caltinay 3670 {
721 caltinay 3691 throw RipleyException("dump() not implemented");
722 caltinay 3670 }
723    
724 caltinay 3691 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
725     {
726     throw RipleyException("borrowSampleReferenceIDs() not implemented");
727     }
728 caltinay 3670
729 caltinay 3691 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
730 caltinay 3670 {
731 caltinay 3691 throw RipleyException("interpolateACross() not implemented");
732 caltinay 3670 }
733    
734 caltinay 3691 bool RipleyDomain::probeInterpolationACross(int fsType_source,
735     const escript::AbstractDomain&, int fsType_target) const
736 caltinay 3670 {
737 caltinay 3691 throw RipleyException("probeInterpolationACross() not implemented");
738 caltinay 3670 }
739    
740 caltinay 3691 void RipleyDomain::setToNormal(escript::Data& normal) const
741 caltinay 3670 {
742 caltinay 3691 throw RipleyException("setToNormal() not implemented");
743 caltinay 3670 }
744    
745 caltinay 3691 void RipleyDomain::setToSize(escript::Data& size) const
746 caltinay 3670 {
747 caltinay 3691 throw RipleyException("setToSize() not implemented");
748 caltinay 3670 }
749    
750 caltinay 3691 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
751 caltinay 3670 {
752 caltinay 3691 throw RipleyException("setToGradient() not implemented");
753 caltinay 3670 }
754    
755 caltinay 3691 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
756 caltinay 3670 {
757 caltinay 3691 throw RipleyException("setToIntegrals() not implemented");
758 caltinay 3670 }
759    
760 caltinay 3691 bool RipleyDomain::ownSample(int fsType, index_t id) const
761 caltinay 3670 {
762 caltinay 3691 throw RipleyException("ownSample() not implemented");
763 caltinay 3670 }
764    
765 caltinay 3691 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
766     const escript::Data& D, const escript::Data& d,
767     const escript::Data& d_dirac, const bool useHRZ) const
768     {
769     throw RipleyException("addPDEToLumpedSystem() not implemented");
770 caltinay 3670 }
771    
772 caltinay 3691 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
773     const escript::Data& Y, const escript::Data& y,
774     const escript::Data& y_contact, const escript::Data& y_dirac) const
775 caltinay 3670 {
776 caltinay 3691 throw RipleyException("addPDEToRHS() not implemented");
777 caltinay 3670 }
778    
779 caltinay 3691 void RipleyDomain::addPDEToTransportProblem(
780     escript::AbstractTransportProblem& tp,
781     escript::Data& source, const escript::Data& M,
782     const escript::Data& A, const escript::Data& B, const escript::Data& C,
783     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
784     const escript::Data& d, const escript::Data& y,
785     const escript::Data& d_contact, const escript::Data& y_contact,
786     const escript::Data& d_dirac, const escript::Data& y_dirac) const
787 caltinay 3670 {
788 caltinay 3691 throw RipleyException("addPDEToTransportProblem() not implemented");
789 caltinay 3670 }
790    
791 caltinay 3691 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
792     const int blocksize, const escript::FunctionSpace& functionspace,
793     const int type) const
794 caltinay 3670 {
795 caltinay 3691 throw RipleyException("newTransportProblem() not implemented");
796 caltinay 3670 }
797    
798 caltinay 3691 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
799     bool reducedColOrder) const
800 caltinay 3670 {
801 caltinay 3691 throw RipleyException("getPattern() not implemented");
802 caltinay 3670 }
803    
804 caltinay 3691 dim_t RipleyDomain::getNumDataPointsGlobal() const
805 caltinay 3670 {
806 caltinay 3691 throw RipleyException("getNumDataPointsGlobal() not implemented");
807 caltinay 3670 }
808    
809 caltinay 3697 IndexVector RipleyDomain::getNumNodesPerDim() const
810     {
811     throw RipleyException("getNumNodesPerDim() not implemented");
812     }
813    
814     IndexVector RipleyDomain::getNumElementsPerDim() const
815     {
816     throw RipleyException("getNumElementsPerDim() not implemented");
817     }
818    
819     IndexVector RipleyDomain::getNumFacesPerBoundary() const
820     {
821     throw RipleyException("getNumFacesPerBoundary() not implemented");
822     }
823    
824     IndexVector RipleyDomain::getNodeDistribution() const
825     {
826     throw RipleyException("getNodeDistribution() not implemented");
827     }
828    
829     pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
830     {
831     throw RipleyException("getFirstCoordAndSpacing() not implemented");
832     }
833    
834 caltinay 3691 dim_t RipleyDomain::getNumFaceElements() const
835 caltinay 3670 {
836 caltinay 3691 throw RipleyException("getNumFaceElements() not implemented");
837 caltinay 3670 }
838    
839 caltinay 3691 dim_t RipleyDomain::getNumElements() const
840 caltinay 3670 {
841 caltinay 3691 throw RipleyException("getNumElements() not implemented");
842 caltinay 3670 }
843    
844 caltinay 3691 dim_t RipleyDomain::getNumNodes() const
845 caltinay 3670 {
846 caltinay 3691 throw RipleyException("getNumNodes() not implemented");
847 caltinay 3670 }
848    
849 caltinay 3691 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
850 caltinay 3670 {
851 caltinay 3691 throw RipleyException("assembleCoordinates() not implemented");
852 caltinay 3670 }
853    
854 caltinay 3748 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
855     const escript::Data& A, const escript::Data& B, const escript::Data& C,
856     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
857     const escript::Data& d, const escript::Data& y) const
858     {
859     throw RipleyException("assemblePDESingle() not implemented");
860     }
861    
862     void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
863     const escript::Data& A, const escript::Data& B, const escript::Data& C,
864     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
865     const escript::Data& d, const escript::Data& y) const
866     {
867     throw RipleyException("assemblePDESystem() not implemented");
868     }
869    
870 caltinay 3711 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
871 caltinay 3701 {
872     throw RipleyException("interpolateNodesOnElements() not implemented");
873     }
874 caltinay 3691
875 caltinay 3711 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
876 caltinay 3701 {
877     throw RipleyException("interpolateNodesOnFaces() not implemented");
878     }
879    
880    
881 caltinay 3670 } // end of namespace ripley
882    

  ViewVC Help
Powered by ViewVC 1.1.26