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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3868 - (show annotations)
Thu Mar 15 06:07:08 2012 UTC (7 years, 1 month ago) by caltinay
File size: 50791 byte(s)
Update to latest trunk

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

  ViewVC Help
Powered by ViewVC 1.1.26