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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3755 - (show annotations)
Thu Jan 5 06:51:31 2012 UTC (7 years, 8 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 32613 byte(s)
Made element IDs globally unique so ownSample can work (to be implemented).
Fixed couple block in 3D so PDEs seem to work now with MPI :-)

1
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 #include <escript/FunctionSpaceFactory.h>
17 #include <pasowrap/SystemMatrixAdapter.h>
18 #include <pasowrap/TransportProblemAdapter.h>
19
20 #include <iomanip>
21
22 using namespace std;
23 using paso::SystemMatrixAdapter;
24 using paso::TransportProblemAdapter;
25
26 namespace ripley {
27
28 escript::Domain_ptr RipleyDomain::loadMesh(const string& filename)
29 {
30 throw RipleyException("loadMesh() not implemented");
31 }
32
33 escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34 {
35 throw RipleyException("readMesh() not implemented");
36 }
37
38 RipleyDomain::RipleyDomain(dim_t dim) :
39 m_numDim(dim),
40 m_status(0)
41 {
42 m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
43 }
44
45 RipleyDomain::~RipleyDomain()
46 {
47 Esys_MPIInfo_free(m_mpiInfo);
48 }
49
50 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 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
62 {
63 switch (fsType) {
64 case DegreesOfFreedom:
65 case ReducedDegreesOfFreedom:
66 case Nodes:
67 case ReducedNodes:
68 case Elements:
69 case ReducedElements:
70 case FaceElements:
71 case ReducedFaceElements:
72 case Points:
73 return true;
74 default:
75 break;
76 }
77 return false;
78 }
79
80 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
81 {
82 switch (fsType) {
83 case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84 case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85 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 break;
94 }
95 return "Invalid function space type code";
96 }
97
98 pair<int,int> RipleyDomain::getDataShape(int fsType) const
99 {
100 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
101 switch (fsType) {
102 case Nodes:
103 case ReducedNodes: //FIXME: reduced
104 return pair<int,int>(1, getNumNodes());
105 case DegreesOfFreedom:
106 case ReducedDegreesOfFreedom: //FIXME: reduced
107 return pair<int,int>(1, getNumDOF());
108 case Elements:
109 return pair<int,int>(ptsPerSample, getNumElements());
110 case FaceElements:
111 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
112 case ReducedElements:
113 return pair<int,int>(1, getNumElements());
114 case ReducedFaceElements:
115 return pair<int,int>(1, getNumFaceElements());
116 case Points:
117 return pair<int,int>(1, 0); //FIXME: dirac
118 default:
119 break;
120 }
121
122 stringstream msg;
123 msg << "getDataShape(): Unsupported function space type " << fsType
124 << " for " << getDescription();
125 throw RipleyException(msg.str());
126 }
127
128 string RipleyDomain::showTagNames() const
129 {
130 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 }
136 return ret.str();
137 }
138
139 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
140 {
141 /*
142 The idea is to use equivalence classes (i.e. types which can be
143 interpolated back and forth):
144 class 0: DOF <-> Nodes
145 class 1: ReducedDOF <-> ReducedNodes
146 class 2: Points
147 class 3: Elements
148 class 4: ReducedElements
149 class 5: FaceElements
150 class 6: ReducedFaceElements
151
152 There is also a set of lines. Interpolation is possible down a line but not
153 between lines.
154 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
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 */
163 if (fs.empty())
164 return false;
165 vector<bool> hasclass(7, false);
166 vector<int> hasline(3, 0);
167 bool hasnodes=false;
168 bool hasrednodes=false;
169 for (size_t i=0; i<fs.size(); ++i) {
170 switch (fs[i]) {
171 case Nodes: hasnodes=true; // fall through
172 case DegreesOfFreedom:
173 hasclass[0]=true;
174 break;
175 case ReducedNodes: hasrednodes=true; // fall through
176 case ReducedDegreesOfFreedom:
177 hasclass[1]=true;
178 break;
179 case Points:
180 hasline[0]=1;
181 hasclass[2]=true;
182 break;
183 case Elements:
184 hasclass[3]=true;
185 hasline[1]=1;
186 break;
187 case ReducedElements:
188 hasclass[4]=true;
189 hasline[1]=1;
190 break;
191 case FaceElements:
192 hasclass[5]=true;
193 hasline[2]=1;
194 break;
195 case ReducedFaceElements:
196 hasclass[6]=true;
197 hasline[2]=1;
198 break;
199 default:
200 return false;
201 }
202 }
203 int numLines=hasline[0]+hasline[1]+hasline[2];
204
205 // fail if we have more than one leaf group
206 // = there are at least two branches we can't interpolate between
207 if (numLines > 1)
208 return false;
209 else if (numLines==1) {
210 // we have points
211 if (hasline[0]==1)
212 resultcode=Points;
213 else if (hasline[1]==1) {
214 if (hasclass[4])
215 resultcode=ReducedElements;
216 else
217 resultcode=Elements;
218 } else { // hasline[2]==1
219 if (hasclass[6])
220 resultcode=ReducedFaceElements;
221 else
222 resultcode=FaceElements;
223 }
224 } else { // numLines==0
225 if (hasclass[1])
226 // something from class 1
227 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228 else // something from class 0
229 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230 }
231 return true;
232 }
233
234 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
235 int fsType_target) const
236 {
237 if (fsType_target != Nodes &&
238 fsType_target != ReducedNodes &&
239 fsType_target != ReducedDegreesOfFreedom &&
240 fsType_target != DegreesOfFreedom &&
241 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 << fsType_target << " for " << getDescription();
249 throw RipleyException(msg.str());
250 }
251
252 switch (fsType_source) {
253 case Nodes:
254 case DegreesOfFreedom:
255 return true;
256 case ReducedNodes:
257 case ReducedDegreesOfFreedom:
258 return (fsType_target != Nodes &&
259 fsType_target != DegreesOfFreedom);
260 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 << fsType_source << " for " << getDescription();
277 throw RipleyException(msg.str());
278 }
279 }
280 }
281
282 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 const int inFS = in.getFunctionSpace().getTypeCode();
299 const int outFS = target.getFunctionSpace().getTypeCode();
300
301 // simplest case: 1:1 copy
302 if (inFS==outFS) {
303 copyData(target, *const_cast<escript::Data*>(&in));
304 // not allowed: reduced->non-reduced
305 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
306 && (outFS==Nodes || outFS==DegreesOfFreedom)) {
307 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 case ReducedNodes: //FIXME: reduced
315 switch (outFS) {
316 case DegreesOfFreedom:
317 case ReducedDegreesOfFreedom: //FIXME: reduced
318 if (getMPISize()==1)
319 copyData(target, *const_cast<escript::Data*>(&in));
320 else
321 nodesToDOF(target,*const_cast<escript::Data*>(&in));
322 break;
323
324 case Nodes:
325 case ReducedNodes: //FIXME: reduced
326 copyData(target, *const_cast<escript::Data*>(&in));
327 break;
328
329 case Elements:
330 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
331 break;
332
333 case ReducedElements:
334 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
335 break;
336
337 case FaceElements:
338 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
339 break;
340
341 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
350 case DegreesOfFreedom:
351 case ReducedDegreesOfFreedom: //FIXME: reduced
352 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 break;
390 default:
391 throw RipleyException(msg.str());
392 }
393 }
394 }
395
396 escript::Data RipleyDomain::getX() const
397 {
398 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 }
428 }
429
430 bool RipleyDomain::isCellOriented(int fsType) const
431 {
432 switch(fsType) {
433 case Nodes:
434 case DegreesOfFreedom:
435 case ReducedDegreesOfFreedom:
436 return false;
437 case Elements:
438 case FaceElements:
439 case Points:
440 case ReducedElements:
441 case ReducedFaceElements:
442 return true;
443 default:
444 break;
445 }
446 stringstream msg;
447 msg << "isCellOriented(): Illegal function space type " << fsType
448 << " on " << getDescription();
449 throw RipleyException(msg.str());
450 }
451
452 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 case DegreesOfFreedom:
462 case ReducedDegreesOfFreedom:
463 case Points:
464 return false;
465 default:
466 break;
467 }
468 stringstream msg;
469 msg << "canTag(): Illegal function space type " << fsType << " on "
470 << getDescription();
471 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 msg << "setTags(): not implemented for "
497 << functionSpaceTypeAsString(fsType);
498 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 msg << "getTagFromSampleNo(): not implemented for "
535 << functionSpaceTypeAsString(fsType);
536 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 msg << "getNumberOfTagsInUse(): not implemented for "
556 << functionSpaceTypeAsString(fsType);
557 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 msg << "borrowListOfTagsInUse(): not implemented for "
576 << functionSpaceTypeAsString(fsType);
577 throw RipleyException(msg.str());
578 }
579 }
580 }
581
582 void RipleyDomain::Print_Mesh_Info(const bool full) const
583 {
584 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
588 // 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 }
597 }
598
599 int RipleyDomain::getSystemMatrixTypeId(const int solver,
600 const int preconditioner, const int package, const bool symmetry) const
601 {
602 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
603 package, symmetry, m_mpiInfo);
604 }
605
606 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
607 const int package, const bool symmetry) const
608 {
609 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
610 package, symmetry, m_mpiInfo);
611 }
612
613 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 {
619 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 throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
625 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
626 if (col_domain!=*this)
627 throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
628 // is the function space type right?
629 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
630 reduceRowOrder=true;
631 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
632 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
633 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
634 reduceColOrder=true;
635 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
636 throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
637
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 }
648
649 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 //TODO: more input param checks (shape, function space etc)
668
669 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
670
671 if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
672 throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
673
674 const int numEq=S->logical_row_block_size;
675 const int numComp=S->logical_col_block_size;
676
677 if (numEq != numComp)
678 throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
679 //TODO: more system matrix checks
680
681 if (numEq==1)
682 assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
683 else
684 assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
685 }
686
687 void RipleyDomain::setNewX(const escript::Data& arg)
688 {
689 throw RipleyException("setNewX(): Operation not supported");
690 }
691
692 //protected
693 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
694 {
695 const dim_t numComp = in.getDataPointSize();
696 const dim_t dpp = in.getNumDataPointsPerSample();
697 out.requireWrite();
698 #pragma omp parallel for
699 for (index_t i=0; i<in.getNumSamples(); i++) {
700 const double* src = in.getSampleDataRO(i);
701 double* dest = out.getSampleDataRW(i);
702 for (index_t c=0; c<numComp; c++) {
703 double res=0.;
704 for (index_t q=0; q<dpp; q++)
705 res+=src[c+q*numComp];
706 *dest++ = res/dpp;
707 }
708 }
709 }
710
711 //protected
712 void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
713 {
714 const dim_t numComp = in.getDataPointSize();
715 out.requireWrite();
716 #pragma omp parallel for
717 for (index_t i=0; i<in.getNumSamples(); i++) {
718 const double* src = in.getSampleDataRO(i);
719 copy(src, src+numComp, out.getSampleDataRW(i));
720 }
721 }
722
723 //protected
724 void RipleyDomain::updateTagsInUse(int fsType) const
725 {
726 IndexVector* tagsInUse=NULL;
727 const IndexVector* tags=NULL;
728 switch(fsType) {
729 case Nodes:
730 tags=&m_nodeTags;
731 tagsInUse=&m_nodeTagsInUse;
732 break;
733 case Elements:
734 case ReducedElements:
735 tags=&m_elementTags;
736 tagsInUse=&m_elementTagsInUse;
737 break;
738 case FaceElements:
739 case ReducedFaceElements:
740 tags=&m_faceTags;
741 tagsInUse=&m_faceTagsInUse;
742 break;
743 default:
744 return;
745 }
746
747 // gather global unique tag values from tags into tagsInUse
748 tagsInUse->clear();
749 index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
750
751 while (true) {
752 // find smallest value bigger than lastFoundValue
753 minFoundValue = INDEX_T_MAX;
754 #pragma omp parallel private(local_minFoundValue)
755 {
756 local_minFoundValue = minFoundValue;
757 #pragma omp for schedule(static) nowait
758 for (size_t i = 0; i < tags->size(); i++) {
759 const index_t v = (*tags)[i];
760 if ((v > lastFoundValue) && (v < local_minFoundValue))
761 local_minFoundValue = v;
762 }
763 #pragma omp critical
764 {
765 if (local_minFoundValue < minFoundValue)
766 minFoundValue = local_minFoundValue;
767 }
768 }
769 #ifdef ESYS_MPI
770 local_minFoundValue = minFoundValue;
771 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
772 #endif
773
774 // if we found a new value add it to the tagsInUse vector
775 if (minFoundValue < INDEX_T_MAX) {
776 tagsInUse->push_back(minFoundValue);
777 lastFoundValue = minFoundValue;
778 } else
779 break;
780 }
781 }
782
783 //
784 // the methods that follow have to be implemented by the subclasses
785 //
786
787 string RipleyDomain::getDescription() const
788 {
789 throw RipleyException("getDescription() not implemented");
790 }
791
792 void RipleyDomain::write(const string& filename) const
793 {
794 throw RipleyException("write() not implemented");
795 }
796
797 void RipleyDomain::dump(const string& filename) const
798 {
799 throw RipleyException("dump() not implemented");
800 }
801
802 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
803 {
804 throw RipleyException("borrowSampleReferenceIDs() not implemented");
805 }
806
807 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
808 {
809 throw RipleyException("interpolateACross() not implemented");
810 }
811
812 bool RipleyDomain::probeInterpolationACross(int fsType_source,
813 const escript::AbstractDomain&, int fsType_target) const
814 {
815 throw RipleyException("probeInterpolationACross() not implemented");
816 }
817
818 void RipleyDomain::setToNormal(escript::Data& normal) const
819 {
820 throw RipleyException("setToNormal() not implemented");
821 }
822
823 void RipleyDomain::setToSize(escript::Data& size) const
824 {
825 throw RipleyException("setToSize() not implemented");
826 }
827
828 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
829 {
830 throw RipleyException("setToGradient() not implemented");
831 }
832
833 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
834 {
835 throw RipleyException("setToIntegrals() not implemented");
836 }
837
838 bool RipleyDomain::ownSample(int fsType, index_t id) const
839 {
840 throw RipleyException("ownSample() not implemented");
841 }
842
843 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
844 const escript::Data& D, const escript::Data& d,
845 const escript::Data& d_dirac, const bool useHRZ) const
846 {
847 throw RipleyException("addPDEToLumpedSystem() not implemented");
848 }
849
850 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
851 const escript::Data& Y, const escript::Data& y,
852 const escript::Data& y_contact, const escript::Data& y_dirac) const
853 {
854 throw RipleyException("addPDEToRHS() not implemented");
855 }
856
857 void RipleyDomain::addPDEToTransportProblem(
858 escript::AbstractTransportProblem& tp,
859 escript::Data& source, const escript::Data& M,
860 const escript::Data& A, const escript::Data& B, const escript::Data& C,
861 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
862 const escript::Data& d, const escript::Data& y,
863 const escript::Data& d_contact, const escript::Data& y_contact,
864 const escript::Data& d_dirac, const escript::Data& y_dirac) const
865 {
866 throw RipleyException("addPDEToTransportProblem() not implemented");
867 }
868
869 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
870 const int blocksize, const escript::FunctionSpace& functionspace,
871 const int type) const
872 {
873 throw RipleyException("newTransportProblem() not implemented");
874 }
875
876 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
877 bool reducedColOrder) const
878 {
879 throw RipleyException("getPattern() not implemented");
880 }
881
882 dim_t RipleyDomain::getNumDataPointsGlobal() const
883 {
884 throw RipleyException("getNumDataPointsGlobal() not implemented");
885 }
886
887 IndexVector RipleyDomain::getNumNodesPerDim() const
888 {
889 throw RipleyException("getNumNodesPerDim() not implemented");
890 }
891
892 IndexVector RipleyDomain::getNumElementsPerDim() const
893 {
894 throw RipleyException("getNumElementsPerDim() not implemented");
895 }
896
897 IndexVector RipleyDomain::getNumFacesPerBoundary() const
898 {
899 throw RipleyException("getNumFacesPerBoundary() not implemented");
900 }
901
902 IndexVector RipleyDomain::getNodeDistribution() const
903 {
904 throw RipleyException("getNodeDistribution() not implemented");
905 }
906
907 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
908 {
909 throw RipleyException("getFirstCoordAndSpacing() not implemented");
910 }
911
912 dim_t RipleyDomain::getNumFaceElements() const
913 {
914 throw RipleyException("getNumFaceElements() not implemented");
915 }
916
917 dim_t RipleyDomain::getNumElements() const
918 {
919 throw RipleyException("getNumElements() not implemented");
920 }
921
922 dim_t RipleyDomain::getNumNodes() const
923 {
924 throw RipleyException("getNumNodes() not implemented");
925 }
926
927 dim_t RipleyDomain::getNumDOF() const
928 {
929 throw RipleyException("getNumDOF() not implemented");
930 }
931
932 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
933 {
934 throw RipleyException("assembleCoordinates() not implemented");
935 }
936
937 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
938 const escript::Data& A, const escript::Data& B, const escript::Data& C,
939 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
940 const escript::Data& d, const escript::Data& y) const
941 {
942 throw RipleyException("assemblePDESingle() not implemented");
943 }
944
945 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
946 const escript::Data& A, const escript::Data& B, const escript::Data& C,
947 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
948 const escript::Data& d, const escript::Data& y) const
949 {
950 throw RipleyException("assemblePDESystem() not implemented");
951 }
952
953 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
954 {
955 throw RipleyException("interpolateNodesOnElements() not implemented");
956 }
957
958 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
959 {
960 throw RipleyException("interpolateNodesOnFaces() not implemented");
961 }
962
963 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
964 {
965 throw RipleyException("nodesToDOF() not implemented");
966 }
967
968 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
969 {
970 throw RipleyException("dofToNodes() not implemented");
971 }
972
973 } // end of namespace ripley
974

  ViewVC Help
Powered by ViewVC 1.1.26