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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4211 - (show annotations)
Mon Feb 18 23:54:46 2013 UTC (6 years, 3 months ago) by caltinay
File size: 52415 byte(s)
Implemented interpolation from Reduced[Face]Elements to [Face]Elements and
changed regularization to compute gradient on Function instead of
ReducedFunction. Results differ slightly so this should help with the accuracy.

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

  ViewVC Help
Powered by ViewVC 1.1.26