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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3773 - (show annotations)
Wed Jan 18 04:27:53 2012 UTC (7 years, 9 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 46774 byte(s)
Adjusted interface for boundary conditions.

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 copyData(target, *const_cast<escript::Data*>(&in));
316 case ReducedNodes: //FIXME: reduced
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);
767 assemblePDEBoundarySingleReduced(S, rhs, escript::Data(),
768 escript::Data(), escript::Data(), d, escript::Data(), y);
769 } else {
770 assemblePDESingle(S, rhs, A, B, C, D, X, Y);
771 assemblePDEBoundarySingle(S, rhs, escript::Data(),
772 escript::Data(), escript::Data(), d, escript::Data(), y);
773 }
774 } else {
775 if (reducedOrder) {
776 assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y);
777 assemblePDEBoundarySystemReduced(S, rhs, escript::Data(),
778 escript::Data(), escript::Data(), d, escript::Data(), y);
779 } else {
780 assemblePDESystem(S, rhs, A, B, C, D, X, Y);
781 assemblePDEBoundarySystem(S, rhs, escript::Data(),
782 escript::Data(), escript::Data(), d, escript::Data(), y);
783 }
784 }
785 }
786
787 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
788 const escript::Data& Y, const escript::Data& y,
789 const escript::Data& y_contact, const escript::Data& y_dirac) const
790 {
791 if (!y_contact.isEmpty())
792 throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
793
794 if (rhs.isEmpty()) {
795 if (!X.isEmpty() || !Y.isEmpty())
796 throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
797 else
798 return;
799 }
800
801 if (rhs.getDataPointSize() == 1) {
802 assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
803 assemblePDEBoundarySingle(NULL, rhs, escript::Data(), y);
804 } else {
805 assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
806 assemblePDEBoundarySystem(NULL, rhs, escript::Data(), y);
807 }
808 }
809
810 void RipleyDomain::setNewX(const escript::Data& arg)
811 {
812 throw RipleyException("setNewX(): Operation not supported");
813 }
814
815 //protected
816 void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
817 {
818 const dim_t numComp = in.getDataPointSize();
819 const dim_t dpp = in.getNumDataPointsPerSample();
820 out.requireWrite();
821 #pragma omp parallel for
822 for (index_t i=0; i<in.getNumSamples(); i++) {
823 const double* src = in.getSampleDataRO(i);
824 double* dest = out.getSampleDataRW(i);
825 for (index_t c=0; c<numComp; c++) {
826 double res=0.;
827 for (index_t q=0; q<dpp; q++)
828 res+=src[c+q*numComp];
829 *dest++ = res/dpp;
830 }
831 }
832 }
833
834 //protected
835 void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
836 {
837 const dim_t numComp = in.getDataPointSize();
838 out.requireWrite();
839 #pragma omp parallel for
840 for (index_t i=0; i<in.getNumSamples(); i++) {
841 const double* src = in.getSampleDataRO(i);
842 copy(src, src+numComp, out.getSampleDataRW(i));
843 }
844 }
845
846 //protected
847 void RipleyDomain::updateTagsInUse(int fsType) const
848 {
849 IndexVector* tagsInUse=NULL;
850 const IndexVector* tags=NULL;
851 switch(fsType) {
852 case Nodes:
853 tags=&m_nodeTags;
854 tagsInUse=&m_nodeTagsInUse;
855 break;
856 case Elements:
857 case ReducedElements:
858 tags=&m_elementTags;
859 tagsInUse=&m_elementTagsInUse;
860 break;
861 case FaceElements:
862 case ReducedFaceElements:
863 tags=&m_faceTags;
864 tagsInUse=&m_faceTagsInUse;
865 break;
866 default:
867 return;
868 }
869
870 // gather global unique tag values from tags into tagsInUse
871 tagsInUse->clear();
872 index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
873
874 while (true) {
875 // find smallest value bigger than lastFoundValue
876 minFoundValue = INDEX_T_MAX;
877 #pragma omp parallel private(local_minFoundValue)
878 {
879 local_minFoundValue = minFoundValue;
880 #pragma omp for schedule(static) nowait
881 for (size_t i = 0; i < tags->size(); i++) {
882 const index_t v = (*tags)[i];
883 if ((v > lastFoundValue) && (v < local_minFoundValue))
884 local_minFoundValue = v;
885 }
886 #pragma omp critical
887 {
888 if (local_minFoundValue < minFoundValue)
889 minFoundValue = local_minFoundValue;
890 }
891 }
892 #ifdef ESYS_MPI
893 local_minFoundValue = minFoundValue;
894 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
895 #endif
896
897 // if we found a new value add it to the tagsInUse vector
898 if (minFoundValue < INDEX_T_MAX) {
899 tagsInUse->push_back(minFoundValue);
900 lastFoundValue = minFoundValue;
901 } else
902 break;
903 }
904 }
905
906 //protected
907 Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
908 const IndexVector& index, const dim_t M, const dim_t N) const
909 {
910 // paso will manage the memory
911 index_t* indexC = MEMALLOC(index.size(), index_t);
912 index_t* ptrC = MEMALLOC(ptr.size(), index_t);
913 copy(index.begin(), index.end(), indexC);
914 copy(ptr.begin(), ptr.end(), ptrC);
915 return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
916 }
917
918 //protected
919 Paso_Pattern* RipleyDomain::createMainPattern() const
920 {
921 IndexVector ptr(1,0);
922 IndexVector index;
923
924 for (index_t i=0; i<getNumDOF(); i++) {
925 // add the DOF itself
926 index.push_back(i);
927 const dim_t num=insertNeighbourNodes(index, i);
928 ptr.push_back(ptr.back()+num+1);
929 }
930
931 return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
932 }
933
934 //protected
935 void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
936 const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
937 {
938 IndexVector ptr(1,0);
939 IndexVector index;
940 for (index_t i=0; i<getNumDOF(); i++) {
941 index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
942 ptr.push_back(ptr.back()+colIndices[i].size());
943 }
944
945 const dim_t M=ptr.size()-1;
946 *colPattern=createPasoPattern(ptr, index, M, N);
947
948 IndexVector rowPtr(1,0);
949 IndexVector rowIndex;
950 for (dim_t id=0; id<N; id++) {
951 dim_t n=0;
952 for (dim_t i=0; i<M; i++) {
953 for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
954 if (index[j]==id) {
955 rowIndex.push_back(i);
956 n++;
957 break;
958 }
959 }
960 }
961 rowPtr.push_back(rowPtr.back()+n);
962 }
963
964 // M and N are now reversed
965 *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
966 }
967
968 //protected
969 void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
970 const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
971 dim_t num_Sol, const vector<double>& array) const
972 {
973 const dim_t numMyCols = mat->pattern->mainPattern->numInput;
974 const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
975 const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
976 const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
977
978 const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
979 const index_t* mainBlock_index = mat->mainBlock->pattern->index;
980 double* mainBlock_val = mat->mainBlock->val;
981 const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
982 const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
983 double* col_coupleBlock_val = mat->col_coupleBlock->val;
984 const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
985 const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
986 double* row_coupleBlock_val = mat->row_coupleBlock->val;
987
988 for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
989 // down columns of array
990 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
991 const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
992 // only look at the matrix rows stored on this processor
993 if (i_row < numMyRows) {
994 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
995 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
996 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
997 if (i_col < numMyCols) {
998 for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
999 if (mainBlock_index[k] == i_col) {
1000 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1001 const dim_t i_Sol=ic+mat->col_block_size*l_col;
1002 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1003 const dim_t i_Eq=ir+mat->row_block_size*l_row;
1004 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())];
1005 }
1006 }
1007 break;
1008 }
1009 }
1010 } else {
1011 for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1012 if (col_coupleBlock_index[k] == i_col - numMyCols) {
1013 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1014 const dim_t i_Sol=ic+mat->col_block_size*l_col;
1015 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1016 const dim_t i_Eq=ir+mat->row_block_size*l_row;
1017 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())];
1018 }
1019 }
1020 break;
1021 }
1022 }
1023 }
1024 }
1025 }
1026 } else {
1027 for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1028 // across rows of array
1029 for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1030 const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1031 if (i_col < numMyCols) {
1032 for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1033 k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1034 {
1035 if (row_coupleBlock_index[k] == i_col) {
1036 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1037 const dim_t i_Sol=ic+mat->col_block_size*l_col;
1038 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1039 const dim_t i_Eq=ir+mat->row_block_size*l_row;
1040 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())];
1041 }
1042 }
1043 break;
1044 }
1045 }
1046 }
1047 }
1048 }
1049 }
1050 }
1051 }
1052 }
1053
1054 //
1055 // the methods that follow have to be implemented by the subclasses
1056 //
1057
1058 string RipleyDomain::getDescription() const
1059 {
1060 throw RipleyException("getDescription() not implemented");
1061 }
1062
1063 void RipleyDomain::write(const string& filename) const
1064 {
1065 throw RipleyException("write() not implemented");
1066 }
1067
1068 void RipleyDomain::dump(const string& filename) const
1069 {
1070 throw RipleyException("dump() not implemented");
1071 }
1072
1073 const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
1074 {
1075 throw RipleyException("borrowSampleReferenceIDs() not implemented");
1076 }
1077
1078 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1079 {
1080 throw RipleyException("interpolateACross() not implemented");
1081 }
1082
1083 bool RipleyDomain::probeInterpolationACross(int fsType_source,
1084 const escript::AbstractDomain&, int fsType_target) const
1085 {
1086 throw RipleyException("probeInterpolationACross() not implemented");
1087 }
1088
1089 void RipleyDomain::setToNormal(escript::Data& normal) const
1090 {
1091 throw RipleyException("setToNormal() not implemented");
1092 }
1093
1094 void RipleyDomain::setToSize(escript::Data& size) const
1095 {
1096 throw RipleyException("setToSize() not implemented");
1097 }
1098
1099 bool RipleyDomain::ownSample(int fsType, index_t id) const
1100 {
1101 throw RipleyException("ownSample() not implemented");
1102 }
1103
1104 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1105 const escript::Data& D, const escript::Data& d,
1106 const escript::Data& d_dirac, const bool useHRZ) const
1107 {
1108 throw RipleyException("addPDEToLumpedSystem() not implemented");
1109 }
1110
1111 void RipleyDomain::addPDEToTransportProblem(
1112 escript::AbstractTransportProblem& tp,
1113 escript::Data& source, const escript::Data& M,
1114 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1115 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1116 const escript::Data& d, const escript::Data& y,
1117 const escript::Data& d_contact, const escript::Data& y_contact,
1118 const escript::Data& d_dirac, const escript::Data& y_dirac) const
1119 {
1120 throw RipleyException("addPDEToTransportProblem() not implemented");
1121 }
1122
1123 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
1124 const int blocksize, const escript::FunctionSpace& functionspace,
1125 const int type) const
1126 {
1127 throw RipleyException("newTransportProblem() not implemented");
1128 }
1129
1130 Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1131 bool reducedColOrder) const
1132 {
1133 throw RipleyException("getPattern() not implemented");
1134 }
1135
1136 dim_t RipleyDomain::getNumDataPointsGlobal() const
1137 {
1138 throw RipleyException("getNumDataPointsGlobal() not implemented");
1139 }
1140
1141 IndexVector RipleyDomain::getNumNodesPerDim() const
1142 {
1143 throw RipleyException("getNumNodesPerDim() not implemented");
1144 }
1145
1146 IndexVector RipleyDomain::getNumElementsPerDim() const
1147 {
1148 throw RipleyException("getNumElementsPerDim() not implemented");
1149 }
1150
1151 IndexVector RipleyDomain::getNumFacesPerBoundary() const
1152 {
1153 throw RipleyException("getNumFacesPerBoundary() not implemented");
1154 }
1155
1156 IndexVector RipleyDomain::getNodeDistribution() const
1157 {
1158 throw RipleyException("getNodeDistribution() not implemented");
1159 }
1160
1161 IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1162 {
1163 throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1164 }
1165
1166 pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1167 {
1168 throw RipleyException("getFirstCoordAndSpacing() not implemented");
1169 }
1170
1171 dim_t RipleyDomain::getNumFaceElements() const
1172 {
1173 throw RipleyException("getNumFaceElements() not implemented");
1174 }
1175
1176 dim_t RipleyDomain::getNumElements() const
1177 {
1178 throw RipleyException("getNumElements() not implemented");
1179 }
1180
1181 dim_t RipleyDomain::getNumNodes() const
1182 {
1183 throw RipleyException("getNumNodes() not implemented");
1184 }
1185
1186 dim_t RipleyDomain::getNumDOF() const
1187 {
1188 throw RipleyException("getNumDOF() not implemented");
1189 }
1190
1191 dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1192 {
1193 throw RipleyException("insertNeighbourNodes() not implemented");
1194 }
1195
1196 void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1197 {
1198 throw RipleyException("assembleCoordinates() not implemented");
1199 }
1200
1201 void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1202 {
1203 throw RipleyException("assembleGradient() not implemented");
1204 }
1205
1206 void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1207 {
1208 throw RipleyException("assembleIntegrate() not implemented");
1209 }
1210
1211 void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1212 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1213 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1214 {
1215 throw RipleyException("assemblePDESingle() not implemented");
1216 }
1217
1218 void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1219 escript::Data& rhs, const escript::Data& a, const escript::Data& b,
1220 const escript::Data& c, const escript::Data& d, const escript::Data& x,
1221 const escript::Data& y) const
1222 {
1223 throw RipleyException("assemblePDEBoundarySingle() not implemented");
1224 }
1225
1226 void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1227 escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1228 const escript::Data& C, const escript::Data& D, const escript::Data& X,
1229 const escript::Data& Y) const
1230 {
1231 throw RipleyException("assemblePDESingleReduced() not implemented");
1232 }
1233
1234 void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1235 escript::Data& rhs, const escript::Data& a, const escript::Data& b,
1236 const escript::Data& c, const escript::Data& d, const escript::Data& x,
1237 const escript::Data& y) const
1238 {
1239 throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1240 }
1241
1242 void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1243 const escript::Data& A, const escript::Data& B, const escript::Data& C,
1244 const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1245 {
1246 throw RipleyException("assemblePDESystem() not implemented");
1247 }
1248
1249 void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1250 escript::Data& rhs, const escript::Data& a, const escript::Data& b,
1251 const escript::Data& c, const escript::Data& d, const escript::Data& x,
1252 const escript::Data& y) const
1253 {
1254 throw RipleyException("assemblePDEBoundarySystem() not implemented");
1255 }
1256
1257 void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1258 escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1259 const escript::Data& C, const escript::Data& D, const escript::Data& X,
1260 const escript::Data& Y) const
1261 {
1262 throw RipleyException("assemblePDESystemReduced() not implemented");
1263 }
1264
1265 void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1266 escript::Data& rhs, const escript::Data& a, const escript::Data& b,
1267 const escript::Data& c, const escript::Data& d, const escript::Data& x,
1268 const escript::Data& y) const
1269 {
1270 throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1271 }
1272
1273 void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1274 {
1275 throw RipleyException("interpolateNodesOnElements() not implemented");
1276 }
1277
1278 void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1279 {
1280 throw RipleyException("interpolateNodesOnFaces() not implemented");
1281 }
1282
1283 void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1284 {
1285 throw RipleyException("nodesToDOF() not implemented");
1286 }
1287
1288 void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1289 {
1290 throw RipleyException("dofToNodes() not implemented");
1291 }
1292
1293 } // end of namespace ripley
1294

  ViewVC Help
Powered by ViewVC 1.1.26