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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3785 - (show annotations)
Wed Jan 25 04:00:33 2012 UTC (7 years, 7 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 45781 byte(s)
Fixed early deallocation of matrix pattern.

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

  ViewVC Help
Powered by ViewVC 1.1.26