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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3776 - (show annotations)
Thu Jan 19 03:48:35 2012 UTC (7 years, 7 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 45966 byte(s)
Finished 2D - all 2D PDE (non-lumped) tests pass now

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

  ViewVC Help
Powered by ViewVC 1.1.26