/[escript]/branches/doubleplusgood/ripley/src/RipleyDomain.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/ripley/src/RipleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4336 - (show annotations)
Thu Mar 21 08:25:18 2013 UTC (6 years, 7 months ago) by jfenwick
File size: 54112 byte(s)
all MEM macros gone except reallocs
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 #include <ripley/RipleyDomain.h>
17 #include <escript/DataFactory.h>
18 #include <escript/FunctionSpaceFactory.h>
19 #include <pasowrap/SystemMatrixAdapter.h>
20 #include <pasowrap/TransportProblemAdapter.h>
21
22 #include <iomanip>
23
24 using namespace std;
25 using paso::SystemMatrixAdapter;
26 using paso::TransportProblemAdapter;
27
28 namespace ripley {
29
30 escript::Domain_ptr RipleyDomain::loadMesh(const string& filename)
31 {
32 throw RipleyException("loadMesh() not implemented");
33 }
34
35 escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
36 {
37 throw RipleyException("readMesh() not implemented");
38 }
39
40 RipleyDomain::RipleyDomain(dim_t dim) :
41 m_numDim(dim),
42 m_status(0)
43 {
44 m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
45 }
46
47 RipleyDomain::~RipleyDomain()
48 {
49 Esys_MPIInfo_free(m_mpiInfo);
50 }
51
52 bool RipleyDomain::operator==(const AbstractDomain& other) const
53 {
54 const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
55 if (o) {
56 return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
57 && m_elementTags==o->m_elementTags
58 && m_faceTags==o->m_faceTags);
59 }
60 return false;
61 }
62
63 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
64 {
65 switch (fsType) {
66 case DegreesOfFreedom:
67 case ReducedDegreesOfFreedom:
68 case Nodes:
69 case ReducedNodes:
70 case Elements:
71 case ReducedElements:
72 case FaceElements:
73 case ReducedFaceElements:
74 case Points:
75 return true;
76 default:
77 break;
78 }
79 return false;
80 }
81
82 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
83 {
84 switch (fsType) {
85 case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
86 case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
87 case Nodes: return "Ripley_Nodes";
88 case ReducedNodes: return "Ripley_ReducedNodes";
89 case Elements: return "Ripley_Elements";
90 case ReducedElements: return "Ripley_ReducedElements";
91 case FaceElements: return "Ripley_FaceElements";
92 case ReducedFaceElements: return "Ripley_ReducedFaceElements";
93 case Points: return "Ripley_Points";
94 default:
95 break;
96 }
97 return "Invalid function space type code";
98 }
99
100 pair<int,int> RipleyDomain::getDataShape(int fsType) const
101 {
102 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
103 switch (fsType) {
104 case Nodes:
105 case ReducedNodes: //FIXME: reduced
106 return pair<int,int>(1, getNumNodes());
107 case DegreesOfFreedom:
108 case ReducedDegreesOfFreedom: //FIXME: reduced
109 return pair<int,int>(1, getNumDOF());
110 case Elements:
111 return pair<int,int>(ptsPerSample, getNumElements());
112 case FaceElements:
113 return pair<int,int>(ptsPerSample/2, getNumFaceElements());
114 case ReducedElements:
115 return pair<int,int>(1, getNumElements());
116 case ReducedFaceElements:
117 return pair<int,int>(1, getNumFaceElements());
118 case Points:
119 return pair<int,int>(1, 0); //FIXME: dirac
120 default:
121 break;
122 }
123
124 stringstream msg;
125 msg << "getDataShape: Invalid function space type " << fsType
126 << " for " << getDescription();
127 throw RipleyException(msg.str());
128 }
129
130 string RipleyDomain::showTagNames() const
131 {
132 stringstream ret;
133 TagMap::const_iterator it;
134 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
135 if (it!=m_tagMap.begin()) ret << ", ";
136 ret << it->first;
137 }
138 return ret.str();
139 }
140
141 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
142 {
143 /*
144 The idea is to use equivalence classes (i.e. types which can be
145 interpolated back and forth):
146 class 0: DOF <-> Nodes
147 class 1: ReducedDOF <-> ReducedNodes
148 class 2: Points
149 class 3: Elements
150 class 4: ReducedElements
151 class 5: FaceElements
152 class 6: ReducedFaceElements
153
154 There is also a set of lines. Interpolation is possible down a line but not
155 between lines.
156 class 0 and 1 belong to all lines so aren't considered.
157 line 0: class 2
158 line 1: class 3,4
159 line 2: class 5,6
160
161 For classes with multiple members (eg class 1) we have vars to record if
162 there is at least one instance. e.g. hasnodes is true if we have at least
163 one instance of Nodes.
164 */
165 if (fs.empty())
166 return false;
167 vector<bool> hasclass(7, false);
168 vector<int> hasline(3, 0);
169 bool hasnodes=false;
170 bool hasrednodes=false;
171 for (size_t i=0; i<fs.size(); ++i) {
172 switch (fs[i]) {
173 case Nodes: hasnodes=true; // fall through
174 case DegreesOfFreedom:
175 hasclass[0]=true;
176 break;
177 case ReducedNodes: hasrednodes=true; // fall through
178 case ReducedDegreesOfFreedom:
179 hasclass[1]=true;
180 break;
181 case Points:
182 hasline[0]=1;
183 hasclass[2]=true;
184 break;
185 case Elements:
186 hasclass[3]=true;
187 hasline[1]=1;
188 break;
189 case ReducedElements:
190 hasclass[4]=true;
191 hasline[1]=1;
192 break;
193 case FaceElements:
194 hasclass[5]=true;
195 hasline[2]=1;
196 break;
197 case ReducedFaceElements:
198 hasclass[6]=true;
199 hasline[2]=1;
200 break;
201 default:
202 return false;
203 }
204 }
205 int numLines=hasline[0]+hasline[1]+hasline[2];
206
207 // fail if we have more than one leaf group
208 // = there are at least two branches we can't interpolate between
209 if (numLines > 1)
210 return false;
211 else if (numLines==1) {
212 // we have points
213 if (hasline[0]==1)
214 resultcode=Points;
215 else if (hasline[1]==1) {
216 if (hasclass[4])
217 resultcode=ReducedElements;
218 else
219 resultcode=Elements;
220 } else { // hasline[2]==1
221 if (hasclass[6])
222 resultcode=ReducedFaceElements;
223 else
224 resultcode=FaceElements;
225 }
226 } else { // numLines==0
227 if (hasclass[1])
228 // something from class 1
229 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
230 else // something from class 0
231 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
232 }
233 return true;
234 }
235
236 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
237 int fsType_target) const
238 {
239 if (!isValidFunctionSpaceType(fsType_target)) {
240 stringstream msg;
241 msg << "probeInterpolationOnDomain: Invalid function space type "
242 << fsType_target << " for " << getDescription();
243 throw RipleyException(msg.str());
244 }
245
246 switch (fsType_source) {
247 case Nodes:
248 case DegreesOfFreedom:
249 return true;
250 case ReducedNodes:
251 case ReducedDegreesOfFreedom:
252 return (fsType_target != Nodes &&
253 fsType_target != DegreesOfFreedom);
254 case Elements:
255 case ReducedElements:
256 return (fsType_target==Elements ||
257 fsType_target==ReducedElements);
258 case FaceElements:
259 case ReducedFaceElements:
260 return (fsType_target==FaceElements ||
261 fsType_target==ReducedFaceElements);
262 case Points:
263 return (fsType_target==fsType_source);
264
265 default: {
266 stringstream msg;
267 msg << "probeInterpolationOnDomain: Invalid function space type "
268 << fsType_source << " for " << getDescription();
269 throw RipleyException(msg.str());
270 }
271 }
272 }
273
274 signed char RipleyDomain::preferredInterpolationOnDomain(int fsType_source,
275 int fsType_target) const
276 {
277 if (!isValidFunctionSpaceType(fsType_target)) {
278 stringstream msg;
279 msg << "preferredInterpolationOnDomain: Invalid function space type "
280 << fsType_target << " for " << getDescription();
281 throw RipleyException(msg.str());
282 }
283
284 if (fsType_source==fsType_target) {
285 return 1;
286 }
287 // There is a complication here in that Nodes and DegreesOfFreedom
288 // can be interpolated to anything, so we need to handle the
289 // reverse case for them specially
290
291 if ((fsType_target==Nodes) || (fsType_target==DegreesOfFreedom))
292 {
293 return -1; // reverse interpolation
294 }
295
296 switch (fsType_source) {
297 case Nodes:
298 case DegreesOfFreedom:
299 return 1;
300 case ReducedNodes:
301 case ReducedDegreesOfFreedom:
302 return (fsType_target != Nodes &&
303 fsType_target != DegreesOfFreedom)?-1:0;
304 case Elements:
305 return (fsType_target==ReducedElements)?1:0;
306 case ReducedElements:
307 return (fsType_target==Elements)?-1:0;
308 case FaceElements:
309 return (fsType_target==ReducedFaceElements)?1:0;
310 case ReducedFaceElements:
311 return (fsType_target==FaceElements)?-1:0;
312 case Points:
313 return false; // other case caught by the if above
314
315 default: {
316 stringstream msg;
317 msg << "probeInterpolationOnDomain: Invalid function space type "
318 << fsType_source << " for " << getDescription();
319 throw RipleyException(msg.str());
320 }
321 }
322 }
323
324 void RipleyDomain::interpolateOnDomain(escript::Data& target,
325 const escript::Data& in) const
326 {
327 const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
328 const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
329 if (inDomain != *this)
330 throw RipleyException("Illegal domain of interpolant");
331 if (targetDomain != *this)
332 throw RipleyException("Illegal domain of interpolation target");
333
334 stringstream msg;
335 msg << "interpolateOnDomain() not implemented for function space "
336 << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
337 << " -> "
338 << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
339
340 const int inFS = in.getFunctionSpace().getTypeCode();
341 const int outFS = target.getFunctionSpace().getTypeCode();
342
343 // simplest case: 1:1 copy
344 if (inFS==outFS) {
345 copyData(target, *const_cast<escript::Data*>(&in));
346 // not allowed: reduced nodes/dof->non-reduced nodes/dof
347 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
348 && (outFS==Nodes || outFS==DegreesOfFreedom)) {
349 throw RipleyException("interpolateOnDomain: Cannot interpolate reduced data to non-reduced data.");
350 } else if ((inFS==Elements && outFS==ReducedElements)
351 || (inFS==FaceElements && outFS==ReducedFaceElements)) {
352 if (in.actsExpanded())
353 averageData(target, *const_cast<escript::Data*>(&in));
354 else
355 copyData(target, *const_cast<escript::Data*>(&in));
356 } else if ((inFS==ReducedElements && outFS==Elements)
357 || (inFS==ReducedFaceElements && outFS==FaceElements)) {
358 multiplyData(target, *const_cast<escript::Data*>(&in));
359 } else {
360 switch (inFS) {
361 case Nodes:
362 case ReducedNodes: //FIXME: reduced
363 switch (outFS) {
364 case DegreesOfFreedom:
365 case ReducedDegreesOfFreedom: //FIXME: reduced
366 if (getMPISize()==1)
367 copyData(target, *const_cast<escript::Data*>(&in));
368 else
369 nodesToDOF(target,*const_cast<escript::Data*>(&in));
370 break;
371
372 case Nodes:
373 case ReducedNodes: //FIXME: reduced
374 copyData(target, *const_cast<escript::Data*>(&in));
375 break;
376
377 case Elements:
378 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
379 break;
380
381 case ReducedElements:
382 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
383 break;
384
385 case FaceElements:
386 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
387 break;
388
389 case ReducedFaceElements:
390 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
391 break;
392
393 default:
394 throw RipleyException(msg.str());
395 }
396 break;
397
398 case DegreesOfFreedom:
399 case ReducedDegreesOfFreedom: //FIXME: reduced
400 switch (outFS) {
401 case Nodes:
402 case ReducedNodes: //FIXME: reduced
403 if (getMPISize()==1)
404 copyData(target, *const_cast<escript::Data*>(&in));
405 else
406 dofToNodes(target, *const_cast<escript::Data*>(&in));
407 break;
408
409 case DegreesOfFreedom:
410 case ReducedDegreesOfFreedom: //FIXME: reduced
411 copyData(target, *const_cast<escript::Data*>(&in));
412 break;
413
414 case Elements:
415 case ReducedElements:
416 if (getMPISize()==1) {
417 interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
418 } else {
419 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
420 escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
421 interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
422 }
423 break;
424
425 case FaceElements:
426 case ReducedFaceElements:
427 if (getMPISize()==1) {
428 interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
429 } else {
430 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
431 escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
432 interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
433 }
434 break;
435
436 default:
437 throw RipleyException(msg.str());
438 }
439 break;
440 default:
441 throw RipleyException(msg.str());
442 }
443 }
444 }
445
446 escript::Data RipleyDomain::getX() const
447 {
448 return escript::continuousFunction(*this).getX();
449 }
450
451 escript::Data RipleyDomain::getNormal() const
452 {
453 return escript::functionOnBoundary(*this).getNormal();
454 }
455
456 escript::Data RipleyDomain::getSize() const
457 {
458 return escript::function(*this).getSize();
459 }
460
461 void RipleyDomain::setToX(escript::Data& arg) const
462 {
463 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
464 *(arg.getFunctionSpace().getDomain()));
465 if (argDomain != *this)
466 throw RipleyException("setToX: Illegal domain of data point locations");
467 if (!arg.isExpanded())
468 throw RipleyException("setToX: Expanded Data object expected");
469
470 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
471 assembleCoordinates(arg);
472 } else {
473 // interpolate the result
474 escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
475 assembleCoordinates(contData);
476 interpolateOnDomain(arg, contData);
477 }
478 }
479
480 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
481 {
482 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
483 *(arg.getFunctionSpace().getDomain()));
484 if (argDomain != *this)
485 throw RipleyException("setToGradient: Illegal domain of gradient argument");
486 const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
487 *(grad.getFunctionSpace().getDomain()));
488 if (gradDomain != *this)
489 throw RipleyException("setToGradient: Illegal domain of gradient");
490
491 switch (grad.getFunctionSpace().getTypeCode()) {
492 case Elements:
493 case ReducedElements:
494 case FaceElements:
495 case ReducedFaceElements:
496 break;
497 default: {
498 stringstream msg;
499 msg << "setToGradient: not supported for "
500 << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
501 throw RipleyException(msg.str());
502 }
503 }
504
505 switch (arg.getFunctionSpace().getTypeCode()) {
506 case DegreesOfFreedom:
507 case ReducedDegreesOfFreedom:
508 case Nodes:
509 case ReducedNodes:
510 break;
511 default: {
512 throw RipleyException("setToGradient: only supported for nodal input data");
513 }
514 }
515
516 if (getMPISize()>1) {
517 if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
518 escript::Data contArg(arg, escript::continuousFunction(*this));
519 assembleGradient(grad, contArg);
520 } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
521 escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
522 assembleGradient(grad, contArg);
523 } else {
524 assembleGradient(grad, *const_cast<escript::Data*>(&arg));
525 }
526 } else {
527 assembleGradient(grad, *const_cast<escript::Data*>(&arg));
528 }
529 }
530
531 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
532 {
533 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
534 *(arg.getFunctionSpace().getDomain()));
535 if (argDomain != *this)
536 throw RipleyException("setToIntegrals: illegal domain of integration kernel");
537
538 switch (arg.getFunctionSpace().getTypeCode()) {
539 case Nodes:
540 case ReducedNodes:
541 case DegreesOfFreedom:
542 case ReducedDegreesOfFreedom:
543 {
544 escript::Data funcArg(arg, escript::function(*this));
545 assembleIntegrate(integrals, funcArg);
546 }
547 break;
548 case Elements:
549 case ReducedElements:
550 case FaceElements:
551 case ReducedFaceElements:
552 assembleIntegrate(integrals, *const_cast<escript::Data*>(&arg));
553 break;
554 default: {
555 stringstream msg;
556 msg << "setToIntegrals: not supported for "
557 << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
558 throw RipleyException(msg.str());
559 }
560 }
561
562 }
563
564 bool RipleyDomain::isCellOriented(int fsType) const
565 {
566 switch(fsType) {
567 case Nodes:
568 case ReducedNodes:
569 case DegreesOfFreedom:
570 case ReducedDegreesOfFreedom:
571 return false;
572 case Elements:
573 case FaceElements:
574 case Points:
575 case ReducedElements:
576 case ReducedFaceElements:
577 return true;
578 default:
579 break;
580 }
581 stringstream msg;
582 msg << "isCellOriented: invalid function space type " << fsType
583 << " on " << getDescription();
584 throw RipleyException(msg.str());
585 }
586
587 bool RipleyDomain::canTag(int fsType) const
588 {
589 switch(fsType) {
590 case Nodes:
591 case Elements:
592 case ReducedElements:
593 case FaceElements:
594 case ReducedFaceElements:
595 return true;
596 case DegreesOfFreedom:
597 case ReducedDegreesOfFreedom:
598 case Points:
599 case ReducedNodes:
600 return false;
601 default:
602 break;
603 }
604 stringstream msg;
605 msg << "canTag: invalid function space type " << fsType << " on "
606 << getDescription();
607 throw RipleyException(msg.str());
608 }
609
610 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
611 {
612 IndexVector* target=NULL;
613 dim_t num=0;
614
615 switch(fsType) {
616 case Nodes:
617 num=getNumNodes();
618 target=&m_nodeTags;
619 break;
620 case Elements:
621 case ReducedElements:
622 num=getNumElements();
623 target=&m_elementTags;
624 break;
625 case FaceElements:
626 case ReducedFaceElements:
627 num=getNumFaceElements();
628 target=&m_faceTags;
629 break;
630 default: {
631 stringstream msg;
632 msg << "setTags: invalid function space type " << fsType;
633 throw RipleyException(msg.str());
634 }
635 }
636 if (target->size() != num) {
637 target->assign(num, -1);
638 }
639
640 escript::Data& mask=*const_cast<escript::Data*>(&cMask);
641 #pragma omp parallel for
642 for (index_t i=0; i<num; i++) {
643 if (mask.getSampleDataRO(i)[0] > 0) {
644 (*target)[i]=newTag;
645 }
646 }
647 updateTagsInUse(fsType);
648 }
649
650 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
651 {
652 switch(fsType) {
653 case Nodes:
654 if (m_nodeTags.size() > sampleNo)
655 return m_nodeTags[sampleNo];
656 break;
657 case Elements:
658 case ReducedElements:
659 if (m_elementTags.size() > sampleNo)
660 return m_elementTags[sampleNo];
661 break;
662 case FaceElements:
663 case ReducedFaceElements:
664 if (m_faceTags.size() > sampleNo)
665 return m_faceTags[sampleNo];
666 break;
667 default: {
668 stringstream msg;
669 msg << "getTagFromSampleNo: invalid function space type " << fsType;
670 throw RipleyException(msg.str());
671 }
672 }
673 return -1;
674 }
675
676 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
677 {
678 switch(fsType) {
679 case Nodes:
680 return m_nodeTagsInUse.size();
681 case Elements:
682 case ReducedElements:
683 return m_elementTagsInUse.size();
684 case FaceElements:
685 case ReducedFaceElements:
686 return m_faceTagsInUse.size();
687 default: {
688 stringstream msg;
689 msg << "getNumberOfTagsInUse: invalid function space type "
690 << fsType;
691 throw RipleyException(msg.str());
692 }
693 }
694 }
695
696 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
697 {
698 switch(fsType) {
699 case Nodes:
700 return &m_nodeTagsInUse[0];
701 case Elements:
702 case ReducedElements:
703 return &m_elementTagsInUse[0];
704 case FaceElements:
705 case ReducedFaceElements:
706 return &m_faceTagsInUse[0];
707 default: {
708 stringstream msg;
709 msg << "borrowListOfTagsInUse: invalid function space type "
710 << fsType;
711 throw RipleyException(msg.str());
712 }
713 }
714 }
715
716 void RipleyDomain::Print_Mesh_Info(const bool full) const
717 {
718 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
719 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
720 cout << "Number of dimensions: " << m_numDim << endl;
721
722 // write tags
723 if (m_tagMap.size() > 0) {
724 cout << "Tags:" << endl;
725 TagMap::const_iterator it;
726 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
727 cout << " " << setw(5) << it->second << " "
728 << it->first << endl;
729 }
730 }
731 }
732
733 int RipleyDomain::getSystemMatrixTypeId(const int solver,
734 const int preconditioner, const int package, const bool symmetry) const
735 {
736 return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
737 package, symmetry, m_mpiInfo);
738 }
739
740 int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
741 const int package, const bool symmetry) const
742 {
743 return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
744 package, symmetry, m_mpiInfo);
745 }
746
747 escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
748 const escript::FunctionSpace& row_functionspace,
749 const int column_blocksize,
750 const escript::FunctionSpace& column_functionspace,
751 const int type) const
752 {
753 bool reduceRowOrder=false;
754 bool reduceColOrder=false;
755 // is the domain right?
756 const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
757 if (row_domain != *this)
758 throw RipleyException("newSystemMatrix: domain of row function space does not match the domain of matrix generator");
759 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
760 if (col_domain != *this)
761 throw RipleyException("newSystemMatrix: domain of column function space does not match the domain of matrix generator");
762 // is the function space type right?
763 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
764 reduceRowOrder=true;
765 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
766 throw RipleyException("newSystemMatrix: illegal function space type for system matrix rows");
767 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
768 reduceColOrder=true;
769 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
770 throw RipleyException("newSystemMatrix: illegal function space type for system matrix columns");
771
772 // generate matrix
773 Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
774 Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
775 row_blocksize, column_blocksize, FALSE);
776 paso::checkPasoError();
777 escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
778 row_functionspace, column_blocksize, column_functionspace));
779 return sma;
780 }
781
782 void RipleyDomain::addPDEToSystem(
783 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
784 const escript::Data& A, const escript::Data& B, const escript::Data& C,
785 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
786 const escript::Data& d, const escript::Data& y,
787 const escript::Data& d_contact, const escript::Data& y_contact,
788 const escript::Data& d_dirac,const escript::Data& y_dirac) const
789 {
790 if (!d_contact.isEmpty() || !y_contact.isEmpty())
791 throw RipleyException("addPDEToSystem: Ripley does not support contact elements");
792
793 paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
794 if (!sma)
795 throw RipleyException("addPDEToSystem: Ripley only accepts Paso system matrices");
796
797 Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
798 assemblePDE(S, rhs, A, B, C, D, X, Y);
799 assemblePDEBoundary(S, rhs, d, y);
800 //assemblePDEDirac(S, rhs, d_dirac, y_dirac);
801
802 }
803
804 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
805 const escript::Data& Y, const escript::Data& y,
806 const escript::Data& y_contact, const escript::Data& y_dirac) const
807 {
808 if (!y_contact.isEmpty())
809 throw RipleyException("addPDEToRHS: Ripley does not support contact elements");
810
811 if (rhs.isEmpty()) {
812 if (!X.isEmpty() || !Y.isEmpty())
813 throw RipleyException("addPDEToRHS: right hand side coefficients are provided but no right hand side vector given");
814 else
815 return;
816 }
817
818 assemblePDE(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
819 assemblePDEBoundary(NULL, rhs, escript::Data(), y);
820 }
821
822 escript::ATP_ptr RipleyDomain::newTransportProblem(const int blocksize,
823 const escript::FunctionSpace& functionspace, const int type) const
824 {
825 bool reduceOrder=false;
826 // is the domain right?
827 const RipleyDomain& domain=dynamic_cast<const RipleyDomain&>(*(functionspace.getDomain()));
828 if (domain != *this)
829 throw RipleyException("newTransportProblem: domain of function space does not match the domain of transport problem generator");
830 // is the function space type right?
831 if (functionspace.getTypeCode()==ReducedDegreesOfFreedom)
832 reduceOrder=true;
833 else if (functionspace.getTypeCode()!=DegreesOfFreedom)
834 throw RipleyException("newTransportProblem: illegal function space type for transport problem");
835
836 // generate matrix
837 Paso_SystemMatrixPattern* pattern=getPattern(reduceOrder, reduceOrder);
838 Paso_TransportProblem* tp = Paso_TransportProblem_alloc(pattern, blocksize);
839 paso::checkPasoError();
840 escript::ATP_ptr atp(new TransportProblemAdapter(tp, blocksize, functionspace));
841 return atp;
842 }
843
844 void RipleyDomain::addPDEToTransportProblem(
845 escript::AbstractTransportProblem& tp,
846 escript::Data& source, const escript::Data& M,
847 const escript::Data& A, const escript::Data& B, const escript::Data& C,
848 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
849 const escript::Data& d, const escript::Data& y,
850 const escript::Data& d_contact, const escript::Data& y_contact,
851 const escript::Data& d_dirac, const escript::Data& y_dirac) const
852 {
853 if (!d_contact.isEmpty() || !y_contact.isEmpty())
854 throw RipleyException("addPDEToTransportProblem: Ripley does not support contact elements");
855
856 paso::TransportProblemAdapter* tpa=dynamic_cast<paso::TransportProblemAdapter*>(&tp);
857 if (!tpa)
858 throw RipleyException("addPDEToTransportProblem: Ripley only accepts Paso transport problems");
859
860 Paso_TransportProblem* ptp = tpa->getPaso_TransportProblem();
861 assemblePDE(ptp->mass_matrix, source, escript::Data(), escript::Data(),
862 escript::Data(), M, escript::Data(), escript::Data());
863 assemblePDE(ptp->transport_matrix, source, A, B, C, D, X, Y);
864 assemblePDEBoundary(ptp->transport_matrix, source, d, y);
865 //assemblePDEDirac(ptp->transport_matrix, source, d_dirac, y_dirac);
866 }
867
868 void RipleyDomain::setNewX(const escript::Data& arg)
869 {
870 throw RipleyException("setNewX(): operation not supported");
871 }
872
873 //protected
874 void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
875 {
876 const dim_t numComp = in.getDataPointSize();
877 out.requireWrite();
878 #pragma omp parallel for
879 for (index_t i=0; i<in.getNumSamples(); i++) {
880 const double* src = in.getSampleDataRO(i);
881 copy(src, src+numComp, out.getSampleDataRW(i));
882 }
883 }
884
885 //protected
886 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
887 {
888 const dim_t numComp = in.getDataPointSize();
889 const dim_t dpp = in.getNumDataPointsPerSample();
890 out.requireWrite();
891 #pragma omp parallel for
892 for (index_t i=0; i<in.getNumSamples(); i++) {
893 const double* src = in.getSampleDataRO(i);
894 double* dest = out.getSampleDataRW(i);
895 for (index_t c=0; c<numComp; c++) {
896 double res=0.;
897 for (index_t q=0; q<dpp; q++)
898 res+=src[c+q*numComp];
899 *dest++ = res/dpp;
900 }
901 }
902 }
903
904 //protected
905 void RipleyDomain::multiplyData(escript::Data& out, escript::Data& in) const
906 {
907 const dim_t numComp = in.getDataPointSize();
908 const dim_t dpp = out.getNumDataPointsPerSample();
909 out.requireWrite();
910 #pragma omp parallel for
911 for (index_t i=0; i<in.getNumSamples(); i++) {
912 const double* src = in.getSampleDataRO(i);
913 double* dest = out.getSampleDataRW(i);
914 for (index_t c=0; c<numComp; c++) {
915 for (index_t q=0; q<dpp; q++)
916 dest[c+q*numComp] = src[c];
917 }
918 }
919 }
920
921 //protected
922 void RipleyDomain::updateTagsInUse(int fsType) const
923 {
924 IndexVector* tagsInUse=NULL;
925 const IndexVector* tags=NULL;
926 switch(fsType) {
927 case Nodes:
928 tags=&m_nodeTags;
929 tagsInUse=&m_nodeTagsInUse;
930 break;
931 case Elements:
932 case ReducedElements:
933 tags=&m_elementTags;
934 tagsInUse=&m_elementTagsInUse;
935 break;
936 case FaceElements:
937 case ReducedFaceElements:
938 tags=&m_faceTags;
939 tagsInUse=&m_faceTagsInUse;
940 break;
941 default:
942 return;
943 }
944
945 // gather global unique tag values from tags into tagsInUse
946 tagsInUse->clear();
947 index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
948
949 while (true) {
950 // find smallest value bigger than lastFoundValue
951 minFoundValue = INDEX_T_MAX;
952 #pragma omp parallel private(local_minFoundValue)
953 {
954 local_minFoundValue = minFoundValue;
955 long i; // should be size_t but omp mutter mutter
956 #pragma omp for schedule(static) private(i) nowait
957 for (i = 0; i < tags->size(); i++) {
958 const index_t v = (*tags)[i];
959 if ((v > lastFoundValue) && (v < local_minFoundValue))
960 local_minFoundValue = v;
961 }
962 #pragma omp critical
963 {
964 if (local_minFoundValue < minFoundValue)
965 minFoundValue = local_minFoundValue;
966 }
967 }
968 #ifdef ESYS_MPI
969 local_minFoundValue = minFoundValue;
970 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
971 #endif
972
973 // if we found a new value add it to the tagsInUse vector
974 if (minFoundValue < INDEX_T_MAX) {
975 tagsInUse->push_back(minFoundValue);
976 lastFoundValue = minFoundValue;
977 } else
978 break;
979 }
980 }
981
982 //protected
983 Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
984 const IndexVector& index, const dim_t M, const dim_t N) const
985 {
986 // paso will manage the memory
987 index_t* indexC = new index_t[index.size()];
988 index_t* ptrC = new index_t[ptr.size()];
989 copy(index.begin(), index.end(), indexC);
990 copy(ptr.begin(), ptr.end(), ptrC);
991 return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
992 }
993
994 //protected
995 Paso_Pattern* RipleyDomain::createMainPattern() const
996 {
997 IndexVector ptr(1,0);
998 IndexVector index;
999
1000 for (index_t i=0; i<getNumDOF(); i++) {
1001 // add the DOF itself
1002 index.push_back(i);
1003 const dim_t num=insertNeighbourNodes(index, i);
1004 ptr.push_back(ptr.back()+num+1);
1005 }
1006
1007 return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
1008 }
1009
1010 //protected
1011 void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
1012 const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
1013 {
1014 IndexVector ptr(1,0);
1015 IndexVector index;
1016 for (index_t i=0; i<getNumDOF(); i++) {
1017 index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
1018 ptr.push_back(ptr.back()+colIndices[i].size());
1019 }
1020
1021 const dim_t M=ptr.size()-1;
1022 *colPattern=createPasoPattern(ptr, index, M, N);
1023
1024 IndexVector rowPtr(1,0);
1025 IndexVector rowIndex;
1026 for (dim_t id=0; id<N; id++) {
1027 dim_t n=0;
1028 for (dim_t i=0; i<M; i++) {
1029 for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
1030 if (index[j]==id) {
1031 rowIndex.push_back(i);
1032 n++;
1033 break;
1034 }
1035 }
1036 }
1037 rowPtr.push_back(rowPtr.back()+n);
1038 }
1039
1040 // M and N are now reversed
1041 *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
1042 }
1043
1044 //protected
1045 void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
1046 const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
1047 dim_t num_Sol, const vector<double>& array) const
1048 {
1049 if (mat->type & MATRIX_FORMAT_TRILINOS_CRS)
1050 throw RipleyException("addToSystemMatrix: TRILINOS_CRS not supported");
1051
1052 const dim_t numMyCols = mat->pattern->mainPattern->numInput;
1053 const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
1054 const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
1055 const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
1056
1057 const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
1058 const index_t* mainBlock_index = mat->mainBlock->pattern->index;
1059 double* mainBlock_val = mat->mainBlock->val;
1060 const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
1061 const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
1062 double* col_coupleBlock_val = mat->col_coupleBlock->val;
1063 const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
1064 const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
1065 double* row_coupleBlock_val = mat->row_coupleBlock->val;
1066 index_t offset=(mat->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
1067
1068 #define UPDATE_BLOCK(VAL) do {\
1069 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {\
1070 const dim_t i_Sol=ic+mat->col_block_size*l_col;\
1071 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {\
1072 const dim_t i_Eq=ir+mat->row_block_size*l_row;\
1073 VAL[k*mat->block_size+ir+mat->row_block_size*ic]\
1074 += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];\
1075 }\
1076 }\
1077 } while(0)
1078
1079 if (mat->type & MATRIX_FORMAT_CSC) {
1080 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1081 // down columns of array
1082 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
1083 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1084 if (i_col < numMyCols) {
1085 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
1086 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
1087 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row+offset;
1088 if (i_row < numMyRows+offset) {
1089 for (dim_t k = mainBlock_ptr[i_col]-offset; k < mainBlock_ptr[i_col+1]-offset; ++k) {
1090 if (mainBlock_index[k] == i_row) {
1091 UPDATE_BLOCK(mainBlock_val);
1092 break;
1093 }
1094 }
1095 } else {
1096 for (dim_t k = col_coupleBlock_ptr[i_col]-offset; k < col_coupleBlock_ptr[i_col+1]-offset; ++k) {
1097 if (row_coupleBlock_index[k] == i_row - numMyRows) {
1098 UPDATE_BLOCK(row_coupleBlock_val);
1099 break;
1100 }
1101 }
1102 }
1103 }
1104 }
1105 } else {
1106 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
1107 // across rows of array
1108 for (dim_t l_row=0; l_row<numSubblocks_Eq; ++l_row) {
1109 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row+offset;
1110 if (i_row < numMyRows+offset) {
1111 for (dim_t k = col_coupleBlock_ptr[i_col-numMyCols]-offset;
1112 k < col_coupleBlock_ptr[i_col-numMyCols+1]-offset; ++k)
1113 {
1114 if (col_coupleBlock_index[k] == i_row) {
1115 UPDATE_BLOCK(col_coupleBlock_val);
1116 break;
1117 }
1118 }
1119 }
1120 }
1121 }
1122 }
1123 }
1124 }
1125 } else {
1126 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
1127 // down columns of array
1128 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
1129 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
1130 // only look at the matrix rows stored on this processor
1131 if (i_row < numMyRows) {
1132 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1133 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
1134 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col+offset;
1135 if (i_col < numMyCols+offset) {
1136 for (dim_t k = mainBlock_ptr[i_row]-offset; k < mainBlock_ptr[i_row+1]-offset; ++k) {
1137 if (mainBlock_index[k] == i_col) {
1138 UPDATE_BLOCK(mainBlock_val);
1139 break;
1140 }
1141 }
1142 } else {
1143 for (dim_t k = col_coupleBlock_ptr[i_row]-offset; k < col_coupleBlock_ptr[i_row+1]-offset; ++k) {
1144 if (col_coupleBlock_index[k] == i_col-numMyCols) {
1145 UPDATE_BLOCK(col_coupleBlock_val);
1146 break;
1147 }
1148 }
1149 }
1150 }
1151 }
1152 } else {
1153 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1154 // across rows of array
1155 for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1156 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col+offset;
1157 if (i_col < numMyCols+offset) {
1158 for (dim_t k = row_coupleBlock_ptr[i_row-numMyRows]-offset;
1159 k < row_coupleBlock_ptr[i_row-numMyRows+1]-offset; ++k)
1160 {
1161 if (row_coupleBlock_index[k] == i_col) {
1162 UPDATE_BLOCK(row_coupleBlock_val);
1163 break;
1164 }
1165 }
1166 }
1167 }
1168 }
1169 }
1170 }
1171 }
1172 }
1173 #undef UPDATE_BLOCK
1174 }
1175
1176 //private
1177 void RipleyDomain::assemblePDE(Paso_SystemMatrix* mat, escript::Data& rhs,
1178 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1179 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1180 {
1181 if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
1182 throw RipleyException("assemblePDE: right hand side coefficients are provided but no right hand side vector given");
1183
1184 vector<int> fsTypes;
1185 if (!A.isEmpty()) fsTypes.push_back(A.getFunctionSpace().getTypeCode());
1186 if (!B.isEmpty()) fsTypes.push_back(B.getFunctionSpace().getTypeCode());
1187 if (!C.isEmpty()) fsTypes.push_back(C.getFunctionSpace().getTypeCode());
1188 if (!D.isEmpty()) fsTypes.push_back(D.getFunctionSpace().getTypeCode());
1189 if (!X.isEmpty()) fsTypes.push_back(X.getFunctionSpace().getTypeCode());
1190 if (!Y.isEmpty()) fsTypes.push_back(Y.getFunctionSpace().getTypeCode());
1191 if (fsTypes.empty())
1192 return;
1193
1194 int fs=fsTypes[0];
1195 if (fs != Elements && fs != ReducedElements)
1196 throw RipleyException("assemblePDE: illegal function space type for coefficients");
1197
1198 for (vector<int>::const_iterator it=fsTypes.begin()+1; it!=fsTypes.end(); it++) {
1199 if (*it != fs) {
1200 throw RipleyException("assemblePDE: coefficient function spaces don't match");
1201 }
1202 }
1203
1204 int numEq, numComp;
1205 if (!mat) {
1206 if (rhs.isEmpty()) {
1207 numEq=numComp=1;
1208 } else {
1209 numEq=numComp=rhs.getDataPointSize();
1210 }
1211 } else {
1212 if (!rhs.isEmpty() && rhs.getDataPointSize()!=mat->logical_row_block_size)
1213 throw RipleyException("assemblePDE: matrix row block size and number of components of right hand side don't match");
1214 numEq = mat->logical_row_block_size;
1215 numComp = mat->logical_col_block_size;
1216 }
1217
1218 if (numEq != numComp)
1219 throw RipleyException("assemblePDE: number of equations and number of solutions don't match");
1220
1221 //TODO: check shape and num samples of coeffs
1222
1223 if (numEq==1) {
1224 if (fs==ReducedElements)
1225 assemblePDESingleReduced(mat, rhs, A, B, C, D, X, Y);
1226 else
1227 assemblePDESingle(mat, rhs, A, B, C, D, X, Y);
1228 } else {
1229 if (fs==ReducedElements)
1230 assemblePDESystemReduced(mat, rhs, A, B, C, D, X, Y);
1231 else
1232 assemblePDESystem(mat, rhs, A, B, C, D, X, Y);
1233 }
1234 }
1235
1236 //private
1237 void RipleyDomain::assemblePDEBoundary(Paso_SystemMatrix* mat,
1238 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1239 {
1240 if (rhs.isEmpty() && !y.isEmpty())
1241 throw RipleyException("assemblePDEBoundary: y provided but no right hand side vector given");
1242
1243 int fs=-1;
1244 if (!d.isEmpty())
1245 fs=d.getFunctionSpace().getTypeCode();
1246 if (!y.isEmpty()) {
1247 if (fs == -1)
1248 fs = y.getFunctionSpace().getTypeCode();
1249 else if (fs != y.getFunctionSpace().getTypeCode())
1250 throw RipleyException("assemblePDEBoundary: coefficient function spaces don't match");
1251 }
1252 if (fs==-1) return;
1253
1254 if (fs != FaceElements && fs != ReducedFaceElements)
1255 throw RipleyException("assemblePDEBoundary: illegal function space type for coefficients");
1256
1257 if (!rhs.isEmpty() && rhs.getDataPointSize() != mat->logical_row_block_size)
1258 throw RipleyException("assemblePDEBoundary: matrix row block size and number of components of right hand side don't match");
1259
1260 const int numEq=mat->logical_row_block_size;
1261 const int numComp=mat->logical_col_block_size;
1262
1263 if (numEq != numComp)
1264 throw RipleyException("assemblePDEBoundary: number of equations and number of solutions don't match");
1265
1266 //TODO: check shape and num samples of coeffs
1267
1268 if (numEq==1) {
1269 if (fs==ReducedFaceElements)
1270 assemblePDEBoundarySingleReduced(mat, rhs, d, y);
1271 else
1272 assemblePDEBoundarySingle(mat, rhs, d, y);
1273 } else {
1274 if (fs==ReducedFaceElements)
1275 assemblePDEBoundarySystemReduced(mat, rhs, d, y);
1276 else
1277 assemblePDEBoundarySystem(mat, rhs, d, y);
1278 }
1279 }
1280
1281 //
1282 // the methods that follow have to be implemented by the subclasses
1283 //
1284
1285 string RipleyDomain::getDescription() const
1286 {
1287 throw RipleyException("getDescription() not implemented");
1288 }
1289
1290 void RipleyDomain::write(const string& filename) const
1291 {
1292 throw RipleyException("write() not implemented");
1293 }
1294
1295 void RipleyDomain::dump(const string& filename) const
1296 {
1297 throw RipleyException("dump() not implemented");
1298 }
1299
1300 void RipleyDomain::readBinaryGrid(escript::Data& out, string filename,
1301 const vector<int>& first, const vector<int>& numValues) const
1302 {
1303 throw RipleyException("readBinaryGrid() not implemented");
1304 }
1305
1306 void RipleyDomain::readNcGrid(escript::Data& out, string filename,
1307 string varname, const vector<int>& first,
1308 const vector<int>& numValues) const
1309 {
1310 throw RipleyException("readNcGrid() not implemented");
1311 }
1312
1313 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
1314 {
1315 throw RipleyException("borrowSampleReferenceIDs() not implemented");
1316 }
1317
1318 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1319 {
1320 throw RipleyException("interpolateACross() not implemented");
1321 }
1322
1323 bool RipleyDomain::probeInterpolationACross(int fsType_source,
1324 const escript::AbstractDomain&, int fsType_target) const
1325 {
1326 throw RipleyException("probeInterpolationACross() not implemented");
1327 }
1328
1329 void RipleyDomain::setToNormal(escript::Data& normal) const
1330 {
1331 throw RipleyException("setToNormal() not implemented");
1332 }
1333
1334 void RipleyDomain::setToSize(escript::Data& size) const
1335 {
1336 throw RipleyException("setToSize() not implemented");
1337 }
1338
1339 bool RipleyDomain::ownSample(int fsType, index_t id) const
1340 {
1341 throw RipleyException("ownSample() not implemented");
1342 }
1343
1344 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1345 bool reducedColOrder) const
1346 {
1347 throw RipleyException("getPattern() not implemented");
1348 }
1349
1350 dim_t RipleyDomain::getNumDataPointsGlobal() const
1351 {
1352 throw RipleyException("getNumDataPointsGlobal() not implemented");
1353 }
1354
1355 IndexVector RipleyDomain::getNumNodesPerDim() const
1356 {
1357 throw RipleyException("getNumNodesPerDim() not implemented");
1358 }
1359
1360 IndexVector RipleyDomain::getNumElementsPerDim() const
1361 {
1362 throw RipleyException("getNumElementsPerDim() not implemented");
1363 }
1364
1365 IndexVector RipleyDomain::getNumFacesPerBoundary() const
1366 {
1367 throw RipleyException("getNumFacesPerBoundary() not implemented");
1368 }
1369
1370 IndexVector RipleyDomain::getNodeDistribution() const
1371 {
1372 throw RipleyException("getNodeDistribution() not implemented");
1373 }
1374
1375 IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1376 {
1377 throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1378 }
1379
1380 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1381 {
1382 throw RipleyException("getFirstCoordAndSpacing() not implemented");
1383 }
1384
1385 dim_t RipleyDomain::getNumFaceElements() const
1386 {
1387 throw RipleyException("getNumFaceElements() not implemented");
1388 }
1389
1390 dim_t RipleyDomain::getNumElements() const
1391 {
1392 throw RipleyException("getNumElements() not implemented");
1393 }
1394
1395 dim_t RipleyDomain::getNumNodes() const
1396 {
1397 throw RipleyException("getNumNodes() not implemented");
1398 }
1399
1400 dim_t RipleyDomain::getNumDOF() const
1401 {
1402 throw RipleyException("getNumDOF() not implemented");
1403 }
1404
1405 dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1406 {
1407 throw RipleyException("insertNeighbourNodes() not implemented");
1408 }
1409
1410 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1411 {
1412 throw RipleyException("assembleCoordinates() not implemented");
1413 }
1414
1415 void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1416 {
1417 throw RipleyException("assembleGradient() not implemented");
1418 }
1419
1420 void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1421 {
1422 throw RipleyException("assembleIntegrate() not implemented");
1423 }
1424
1425 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1426 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1427 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1428 {
1429 throw RipleyException("assemblePDESingle() not implemented");
1430 }
1431
1432 void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1433 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1434 {
1435 throw RipleyException("assemblePDEBoundarySingle() not implemented");
1436 }
1437
1438 void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1439 escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1440 const escript::Data& C, const escript::Data& D, const escript::Data& X,
1441 const escript::Data& Y) const
1442 {
1443 throw RipleyException("assemblePDESingleReduced() not implemented");
1444 }
1445
1446 void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1447 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1448 {
1449 throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1450 }
1451
1452 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1453 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1454 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1455 {
1456 throw RipleyException("assemblePDESystem() not implemented");
1457 }
1458
1459 void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1460 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1461 {
1462 throw RipleyException("assemblePDEBoundarySystem() not implemented");
1463 }
1464
1465 void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1466 escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1467 const escript::Data& C, const escript::Data& D, const escript::Data& X,
1468 const escript::Data& Y) const
1469 {
1470 throw RipleyException("assemblePDESystemReduced() not implemented");
1471 }
1472
1473 void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1474 escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1475 {
1476 throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1477 }
1478
1479 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1480 {
1481 throw RipleyException("interpolateNodesOnElements() not implemented");
1482 }
1483
1484 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1485 {
1486 throw RipleyException("interpolateNodesOnFaces() not implemented");
1487 }
1488
1489 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1490 {
1491 throw RipleyException("nodesToDOF() not implemented");
1492 }
1493
1494 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1495 {
1496 throw RipleyException("dofToNodes() not implemented");
1497 }
1498
1499 } // end of namespace ripley
1500

  ViewVC Help
Powered by ViewVC 1.1.26