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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3764 - (show annotations)
Tue Jan 10 06:31:15 2012 UTC (7 years, 8 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 44747 byte(s)
Fixed integration now that elements overlap + code cleanup.

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

  ViewVC Help
Powered by ViewVC 1.1.26