/[escript]/trunk/speckley/src/SpeckleyDomain.cpp
ViewVC logotype

Contents of /trunk/speckley/src/SpeckleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5478 - (show annotations)
Wed Feb 18 01:57:49 2015 UTC (4 years, 7 months ago) by sshaw
File size: 34874 byte(s)
introducing shiny new reduced function spaces to speckley
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2015 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 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #define ESNEEDPYTHON
18 #include "esysUtils/first.h"
19
20
21 #include <speckley/SpeckleyDomain.h>
22 #include <escript/DataFactory.h>
23 #include <escript/FunctionSpaceFactory.h>
24 #include <speckley/domainhelpers.h>
25
26 #include <iomanip>
27
28 namespace bp = boost::python;
29
30 using namespace std;
31
32 namespace speckley {
33
34 void tupleListToMap(DataMap& mapping, const bp::list& list)
35 {
36 for (int i = 0; i < len(list); i++) {
37 if (!bp::extract<bp::tuple>(list[i]).check())
38 throw SpeckleyException("Passed in list contains objects"
39 " other than tuples");
40 bp::tuple t = bp::extract<bp::tuple>(list[i]);
41 if (len(t) != 2 || !bp::extract<string>(t[0]).check() ||
42 !bp::extract<escript::Data>(t[1]).check())
43 throw SpeckleyException("The passed in list must contain tuples"
44 " of the form (string, escript::data)");
45 mapping[bp::extract<string>(t[0])] = bp::extract<escript::Data>(t[1]);
46 }
47 }
48
49 SpeckleyDomain::SpeckleyDomain(dim_t dim, int order, escript::SubWorld_ptr p) :
50 m_numDim(dim),
51 m_status(0),
52 m_order(order)
53 {
54 if (p.get() == NULL)
55 m_mpiInfo = esysUtils::makeInfo(MPI_COMM_WORLD);
56 else
57 m_mpiInfo = p->getMPI();
58
59 assembler_type = DEFAULT_ASSEMBLER;
60 }
61
62 SpeckleyDomain::~SpeckleyDomain()
63 {
64 // cleanup of MPI is dealt with by shared_ptr
65 }
66
67 bool SpeckleyDomain::operator==(const AbstractDomain& other) const
68 {
69 const SpeckleyDomain* o=dynamic_cast<const SpeckleyDomain*>(&other);
70 if (o) {
71 return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
72 && m_elementTags==o->m_elementTags);
73 }
74 return false;
75 }
76
77 bool SpeckleyDomain::isValidFunctionSpaceType(int fsType) const
78 {
79 switch (fsType) {
80 case DegreesOfFreedom:
81 case Nodes:
82 case Elements:
83 case ReducedElements:
84 case Points:
85 return true;
86 default:
87 break;
88 }
89 return false;
90 }
91
92 string SpeckleyDomain::functionSpaceTypeAsString(int fsType) const
93 {
94 switch (fsType) {
95 case DegreesOfFreedom: return "Speckley_DegreesOfFreedom [Solution(domain)]";
96 case ReducedDegreesOfFreedom: return "Speckley_ReducedDegreesOfFreedom [ReducedSolution(domain)]";
97 case Nodes: return "Speckley_Nodes [ContinuousFunction(domain)]";
98 case ReducedNodes: return "Speckley_ReducedNodes [ReducedContinuousFunction(domain)]";
99 case Elements: return "Speckley_Elements [Function(domain)]";
100 case ReducedElements: return "Speckley_ReducedElements [ReducedFunction(domain)]";
101 case FaceElements: return "Speckley_FaceElements [FunctionOnBoundary(domain)]";
102 case ReducedFaceElements: return "Speckley_ReducedFaceElements [ReducedFunctionOnBoundary(domain)]";
103 case Points: return "Speckley_Points [DiracDeltaFunctions(domain)]";
104 default:
105 break;
106 }
107 return "Invalid function space type code";
108 }
109
110 pair<int,dim_t> SpeckleyDomain::getDataShape(int fsType) const
111 {
112 int ptsPerSample = (m_order + 1) * (m_order + 1);
113 if (m_numDim==3)
114 ptsPerSample *= (m_order + 1);
115 switch (fsType) {
116 case Nodes:
117 return pair<int,dim_t>(1, getNumNodes());
118 case DegreesOfFreedom:
119 return pair<int,dim_t>(1, getNumDOF());
120 case Elements:
121 return pair<int,dim_t>(ptsPerSample, getNumElements());
122 case ReducedElements:
123 return pair<int,dim_t>(1, getNumElements());
124 case Points:
125 return pair<int,dim_t>(1, m_diracPoints.size());
126 default:
127 break;
128 }
129
130 stringstream msg;
131 msg << "getDataShape: Invalid function space type " << fsType
132 << " for " << getDescription();
133 throw SpeckleyException(msg.str());
134 }
135
136 string SpeckleyDomain::showTagNames() const
137 {
138 stringstream ret;
139 TagMap::const_iterator it;
140 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
141 if (it!=m_tagMap.begin()) ret << ", ";
142 ret << it->first;
143 }
144 return ret.str();
145 }
146
147 bool SpeckleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
148 {
149 /*
150 The idea is to use equivalence classes (i.e. types which can be
151 interpolated back and forth):
152 class 0: DOF <-> Nodes
153 class 1: ReducedDOF <-> ReducedNodes
154 class 2: Points
155 class 3: Elements
156 class 4: ReducedElements
157
158 There is also a set of lines. Interpolation is possible down a line but not
159 between lines.
160 class 0 and 1 belong to all lines so aren't considered.
161 line 0: class 2
162 line 1: class 3,4
163
164 For classes with multiple members (eg class 1) we have vars to record if
165 there is at least one instance. e.g. hasnodes is true if we have at least
166 one instance of Nodes.
167 */
168 if (fs.empty())
169 return false;
170 vector<bool> hasclass(7, false);
171 vector<int> hasline(3, 0);
172 bool hasnodes=false;
173 bool hasrednodes=false;
174 for (size_t i=0; i<fs.size(); ++i) {
175 switch (fs[i]) {
176 case Nodes: hasnodes=true; // fall through
177 case DegreesOfFreedom:
178 hasclass[0]=true;
179 break;
180 case ReducedNodes: hasrednodes=true; // fall through
181 case ReducedDegreesOfFreedom:
182 hasclass[1]=true;
183 break;
184 case Points:
185 hasline[0]=1;
186 hasclass[2]=true;
187 break;
188 case Elements:
189 hasclass[3]=true;
190 hasline[1]=1;
191 break;
192 case ReducedElements:
193 hasclass[4]=true;
194 hasline[1]=1;
195 break;
196 default:
197 return false;
198 }
199 }
200 int numLines=hasline[0]+hasline[1];
201
202 // fail if we have more than one leaf group
203 // = there are at least two branches we can't interpolate between
204 if (numLines > 1)
205 return false;
206 else if (numLines==1) {
207 // we have points
208 if (hasline[0]==1)
209 resultcode=Points;
210 else if (hasline[1]==1) {
211 if (hasclass[4])
212 resultcode=ReducedElements;
213 else
214 resultcode=Elements;
215 }
216 } else { // numLines==0
217 if (hasclass[1])
218 // something from class 1
219 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
220 else // something from class 0
221 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
222 }
223 return true;
224 }
225
226 bool SpeckleyDomain::probeInterpolationOnDomain(int fsType_source,
227 int fsType_target) const
228 {
229 if (!isValidFunctionSpaceType(fsType_target)) {
230 stringstream msg;
231 msg << "probeInterpolationOnDomain: Invalid function space type "
232 << fsType_target << " for " << getDescription();
233 throw SpeckleyException(msg.str());
234 }
235
236 switch (fsType_source) {
237 case Nodes:
238 case DegreesOfFreedom:
239 return true;
240 case ReducedNodes:
241 case ReducedDegreesOfFreedom:
242 return (fsType_target != Nodes &&
243 fsType_target != DegreesOfFreedom);
244 case Elements:
245 return (fsType_target==Elements || fsType_target==ReducedElements || fsType_target==Nodes);
246 case Points:
247 return (fsType_target==fsType_source);
248 case ReducedElements:
249 return (fsType_target==Elements || fsType_target==Nodes);
250
251 default: {
252 stringstream msg;
253 msg << "probeInterpolationOnDomain: Invalid function space type "
254 << fsType_source << " for " << getDescription();
255 throw SpeckleyException(msg.str());
256 }
257 }
258 }
259
260 signed char SpeckleyDomain::preferredInterpolationOnDomain(int fsType_source,
261 int fsType_target) const
262 {
263 if (!isValidFunctionSpaceType(fsType_target)) {
264 stringstream msg;
265 msg << "preferredInterpolationOnDomain: Invalid function space type "
266 << fsType_target << " for " << getDescription();
267 throw SpeckleyException(msg.str());
268 }
269
270 if (fsType_source==fsType_target) {
271 return 1;
272 }
273 // There is a complication here in that Nodes and DegreesOfFreedom
274 // can be interpolated to anything, so we need to handle the
275 // reverse case for them specially
276
277 if ((fsType_target==Nodes) || (fsType_target==DegreesOfFreedom)) {
278 return -1; // reverse interpolation
279 }
280
281 switch (fsType_source) {
282 case Nodes:
283 case DegreesOfFreedom:
284 return 1;
285 case ReducedNodes:
286 case ReducedDegreesOfFreedom:
287 return (fsType_target != Nodes &&
288 fsType_target != DegreesOfFreedom) ? -1 : 0;
289 case Elements:
290 return (fsType_target==ReducedElements) ? -1 : 0;
291 case ReducedElements:
292 return (fsType_target==Elements) ? 1 : 0;
293 case Points:
294 return 0; // other case caught by the if above
295 default: {
296 stringstream msg;
297 msg << "probeInterpolationOnDomain: Invalid function space type "
298 << fsType_source << " for " << getDescription();
299 throw SpeckleyException(msg.str());
300 }
301 }
302 }
303
304 void SpeckleyDomain::interpolateOnDomain(escript::Data& target,
305 const escript::Data& in) const
306 {
307 const SpeckleyDomain& inDomain=dynamic_cast<const SpeckleyDomain&>(*(in.getFunctionSpace().getDomain()));
308 const SpeckleyDomain& targetDomain=dynamic_cast<const SpeckleyDomain&>(*(target.getFunctionSpace().getDomain()));
309 if (inDomain != *this)
310 throw SpeckleyException("Illegal domain of interpolant");
311 if (targetDomain != *this)
312 throw SpeckleyException("Illegal domain of interpolation target");
313
314 stringstream msg;
315 msg << "interpolateOnDomain() not implemented for function space "
316 << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
317 << " -> "
318 << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
319
320 const int inFS = in.getFunctionSpace().getTypeCode();
321 const int outFS = target.getFunctionSpace().getTypeCode();
322
323 // simplest case: 1:1 copy
324 if (inFS==outFS) {
325 copyData(target, in);
326 } else if ((inFS==Elements || inFS==ReducedElements) && outFS==Nodes) {
327 interpolateElementsOnNodes(target, in);
328 } else if (inFS==ReducedElements && outFS==Elements) {
329 multiplyData(target, in);
330 } else if (inFS==Elements && outFS==ReducedElements) {
331 reduceElements(target, in);
332 } else {
333 switch (inFS) {
334 case Nodes:
335 switch (outFS) {
336 case DegreesOfFreedom:
337 case Nodes:
338 copyData(target, in);
339 break;
340 case Elements:
341 interpolateNodesOnElements(target, in, false);
342 break;
343 case ReducedElements:
344 interpolateNodesOnElements(target, in, true);
345 break;
346 case Points:
347 {
348 const dim_t numComp = in.getDataPointSize();
349 const int nDirac = m_diracPoints.size();
350 target.requireWrite();
351 #pragma omp parallel for
352 for (int i = 0; i < nDirac; i++) {
353 const double* src = in.getSampleDataRO(m_diracPoints[i].node);
354 copy(src, src+numComp, target.getSampleDataRW(i));
355 }
356 }
357 break;
358 default:
359 throw SpeckleyException(msg.str());
360 }
361 break;
362
363 case DegreesOfFreedom:
364 switch (outFS) {
365 case Nodes:
366 case DegreesOfFreedom:
367 copyData(target, in);
368 break;
369
370 case Elements:
371 case ReducedElements:
372 if (getMPISize()==1) {
373 interpolateNodesOnElements(target, in, outFS==ReducedElements);
374 } else {
375 escript::Data contIn(in, escript::continuousFunction(*this));
376 interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
377 }
378 break;
379
380 default:
381 throw SpeckleyException(msg.str());
382 }
383 break;
384 case Points:
385 switch(outFS) {
386 case Nodes:
387 {
388 const dim_t numComp = in.getDataPointSize();
389 const int nDirac = m_diracPoints.size();
390 target.requireWrite();
391 #pragma omp parallel for
392 for (int i = 0; i < nDirac; i++) {
393 const double* src = in.getSampleDataRO(i);
394 copy(src, src+numComp, target.getSampleDataRW(m_diracPoints[i].node));
395 }
396 }
397
398 }
399 break;
400 default:
401 throw SpeckleyException(msg.str());
402 }
403 }
404 }
405
406 escript::Data SpeckleyDomain::getX() const
407 {
408 return escript::continuousFunction(*this).getX();
409 }
410
411 escript::Data SpeckleyDomain::getNormal() const
412 {
413 throw SpeckleyException("Speckley doesn't support getNormal");
414 }
415
416 escript::Data SpeckleyDomain::getSize() const
417 {
418 return escript::function(*this).getSize();
419 }
420
421 void SpeckleyDomain::setToX(escript::Data& arg) const
422 {
423 const SpeckleyDomain& argDomain=dynamic_cast<const SpeckleyDomain&>(
424 *(arg.getFunctionSpace().getDomain()));
425 if (argDomain != *this)
426 throw SpeckleyException("setToX: Illegal domain of data point locations");
427 if (!arg.isExpanded())
428 throw SpeckleyException("setToX: Expanded Data object expected");
429
430 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
431 assembleCoordinates(arg);
432 } else {
433 // interpolate the result
434 escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
435 assembleCoordinates(contData);
436 interpolateOnDomain(arg, contData);
437 }
438 }
439
440 void SpeckleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
441 {
442 const SpeckleyDomain& argDomain=dynamic_cast<const SpeckleyDomain&>(
443 *(arg.getFunctionSpace().getDomain()));
444 if (argDomain != *this)
445 throw SpeckleyException("setToGradient: Illegal domain of gradient argument");
446 const SpeckleyDomain& gradDomain=dynamic_cast<const SpeckleyDomain&>(
447 *(grad.getFunctionSpace().getDomain()));
448 if (gradDomain != *this)
449 throw SpeckleyException("setToGradient: Illegal domain of gradient");
450
451 switch (grad.getFunctionSpace().getTypeCode()) {
452 case Nodes:
453 case Elements:
454 case ReducedElements:
455 break;
456 default: {
457 stringstream msg;
458 msg << "setToGradient: not supported for "
459 << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
460 throw SpeckleyException(msg.str());
461 }
462 }
463
464 switch (arg.getFunctionSpace().getTypeCode()) {
465 case DegreesOfFreedom:
466 case Nodes:
467 case Elements:
468 break;
469 default: {
470 throw SpeckleyException("setToGradient: only supported for nodal input data");
471 }
472 }
473
474 if (getMPISize() > 1) {
475 if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
476 escript::Data contArg(arg, escript::continuousFunction(*this));
477 assembleGradient(grad, contArg);
478 } else {
479 assembleGradient(grad, arg);
480 }
481 } else {
482 assembleGradient(grad, arg);
483 }
484 }
485
486 void SpeckleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
487 {
488 const SpeckleyDomain& argDomain=dynamic_cast<const SpeckleyDomain&>(
489 *(arg.getFunctionSpace().getDomain()));
490 if (argDomain != *this)
491 throw SpeckleyException("setToIntegrals: illegal domain of integration kernel");
492
493 switch (arg.getFunctionSpace().getTypeCode()) {
494 case Nodes:
495 case DegreesOfFreedom:
496 {
497 escript::Data funcArg(arg, escript::function(*this));
498 assembleIntegrate(integrals, funcArg);
499 }
500 break;
501 case Elements:
502 case ReducedElements:
503 assembleIntegrate(integrals, arg);
504 break;
505 default: {
506 stringstream msg;
507 msg << "setToIntegrals: not supported for "
508 << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
509 throw SpeckleyException(msg.str());
510 }
511 }
512
513 }
514
515 bool SpeckleyDomain::isCellOriented(int fsType) const
516 {
517 switch(fsType) {
518 case Nodes:
519 case DegreesOfFreedom:
520 return false;
521 case Elements:
522 case ReducedElements:
523 case Points:
524 return true;
525 default:
526 break;
527 }
528 stringstream msg;
529 msg << "isCellOriented: invalid function space type " << fsType
530 << " on " << getDescription();
531 throw SpeckleyException(msg.str());
532 }
533
534 bool SpeckleyDomain::canTag(int fsType) const
535 {
536 switch(fsType) {
537 case Nodes:
538 case Elements:
539 case ReducedElements:
540 case Points:
541 return true;
542 case DegreesOfFreedom:
543 return false;
544 default:
545 break;
546 }
547 stringstream msg;
548 msg << "canTag: invalid function space type " << fsType << " on "
549 << getDescription();
550 throw SpeckleyException(msg.str());
551 }
552
553 void SpeckleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const
554 {
555 vector<int>* target=NULL;
556 dim_t num=0;
557
558 switch(fsType) {
559 case Nodes:
560 num=getNumNodes();
561 target=&m_nodeTags;
562 break;
563 case Elements:
564 num=getNumElements();
565 target=&m_elementTags;
566 break;
567 default: {
568 stringstream msg;
569 msg << "setTags: invalid function space type " << fsType;
570 throw SpeckleyException(msg.str());
571 }
572 }
573 if (target->size() != num) {
574 target->assign(num, -1);
575 }
576
577 #pragma omp parallel for
578 for (index_t i=0; i<num; i++) {
579 if (mask.getSampleDataRO(i)[0] > 0) {
580 (*target)[i]=newTag;
581 }
582 }
583 updateTagsInUse(fsType);
584 }
585
586 int SpeckleyDomain::getTagFromSampleNo(int fsType, index_t sampleNo) const
587 {
588 switch(fsType) {
589 case Nodes:
590 if (m_nodeTags.size() > sampleNo)
591 return m_nodeTags[sampleNo];
592 break;
593 case Elements:
594 case ReducedElements:
595 if (m_elementTags.size() > sampleNo)
596 return m_elementTags[sampleNo];
597 break;
598 case Points:
599 if (m_diracPoints.size() > sampleNo)
600 return m_diracPoints[sampleNo].tag;
601 break;
602 default: {
603 stringstream msg;
604 msg << "getTagFromSampleNo: invalid function space type " << fsType;
605 throw SpeckleyException(msg.str());
606 }
607 }
608 return -1;
609 }
610
611 int SpeckleyDomain::getNumberOfTagsInUse(int fsType) const
612 {
613 switch(fsType) {
614 case Nodes:
615 return m_nodeTagsInUse.size();
616 case Elements:
617 case ReducedElements:
618 return m_elementTagsInUse.size();
619 default: {
620 stringstream msg;
621 msg << "getNumberOfTagsInUse: invalid function space type "
622 << fsType;
623 throw SpeckleyException(msg.str());
624 }
625 }
626 }
627
628 const int* SpeckleyDomain::borrowListOfTagsInUse(int fsType) const
629 {
630 switch(fsType) {
631 case Nodes:
632 return &m_nodeTagsInUse[0];
633 case Elements:
634 case ReducedElements:
635 return &m_elementTagsInUse[0];
636 default: {
637 stringstream msg;
638 msg << "borrowListOfTagsInUse: invalid function space type "
639 << fsType;
640 throw SpeckleyException(msg.str());
641 }
642 }
643 }
644
645 void SpeckleyDomain::Print_Mesh_Info(bool full) const
646 {
647 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
648 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
649 cout << "Number of dimensions: " << m_numDim << endl;
650 cout << "Number of elements per rank: " << getNumElements() << endl;
651 // write tags
652 if (m_tagMap.size() > 0) {
653 cout << "Tags:" << endl;
654 TagMap::const_iterator it;
655 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
656 cout << " " << setw(5) << it->second << " " << it->first << endl;
657 }
658 }
659 }
660
661 int SpeckleyDomain::getSystemMatrixTypeId(const boost::python::object& options) const
662 {
663 throw SpeckleyException("System matrices not supported by Speckley");
664 }
665
666 int SpeckleyDomain::getTransportTypeId(int solver, int preconditioner,
667 int package, bool symmetry) const
668 {
669 throw SpeckleyException("Transport problems not supported by Speckley");
670 }
671
672 escript::ASM_ptr SpeckleyDomain::newSystemMatrix(int row_blocksize,
673 const escript::FunctionSpace& row_functionspace, int column_blocksize,
674 const escript::FunctionSpace& column_functionspace, int type) const
675 {
676 throw SpeckleyException("Speckley domains do not support system matrices");
677 }
678
679 void SpeckleyDomain::addToSystem(escript::AbstractSystemMatrix& mat,
680 escript::Data& rhs, const DataMap& coefs,
681 Assembler_ptr assembler) const
682 {
683 throw SpeckleyException("Speckley domains do not support system matrices");
684 }
685
686 void SpeckleyDomain::addToSystemFromPython(escript::AbstractSystemMatrix& mat,
687 escript::Data& rhs,
688 const bp::list& data,
689 Assembler_ptr assembler) const
690 {
691 DataMap mapping;
692 tupleListToMap(mapping, data);
693 addToSystem(mat, rhs, mapping, assembler);
694 }
695
696 Assembler_ptr SpeckleyDomain::createAssemblerFromPython(const string type,
697 const bp::list& options) const
698 {
699 DataMap mapping;
700 tupleListToMap(mapping, options);
701 return createAssembler(type, mapping);
702 }
703
704 void SpeckleyDomain::addToRHSFromPython(escript::Data& rhs,
705 const bp::list& data, Assembler_ptr assembler) const
706 {
707 DataMap mapping;
708 tupleListToMap(mapping, data);
709 rhs.expand(); //yes, this is absolutely neccessary
710 addToRHS(rhs, mapping, assembler);
711 }
712
713 void SpeckleyDomain::addToRHS(escript::Data& rhs, const DataMap& coefs,
714 Assembler_ptr assembler) const
715 {
716 if (isNotEmpty("y_contact", coefs))
717 throw SpeckleyException(
718 "addPDEToRHS: Speckley does not support contact elements");
719
720 if (rhs.isEmpty()) {
721 if (isNotEmpty("X", coefs) || isNotEmpty("Y", coefs))
722 throw SpeckleyException(
723 "addPDEToRHS: right hand side coefficients are provided "
724 "but no right hand side vector given");
725 else
726 return;
727 }
728
729 assemblePDE(NULL, rhs, coefs, assembler);
730 assemblePDEBoundary(NULL, rhs, coefs, assembler);
731 assemblePDEDirac(NULL, rhs, coefs, assembler);
732 }
733
734 escript::ATP_ptr SpeckleyDomain::newTransportProblem(int blocksize,
735 const escript::FunctionSpace& functionspace, int type) const
736 {
737 throw SpeckleyException("Speckley domains do not support transport problems");
738 }
739
740 void SpeckleyDomain::addPDEToTransportProblemFromPython(
741 escript::AbstractTransportProblem& tp, escript::Data& source,
742 const bp::list& data, Assembler_ptr assembler) const
743 {
744 throw SpeckleyException("Speckley domains do not support transport problems");
745 }
746
747 void SpeckleyDomain::addPDEToTransportProblem(
748 escript::AbstractTransportProblem& tp,
749 escript::Data& source, const DataMap& coefs,
750 Assembler_ptr assembler) const
751 {
752 throw SpeckleyException("Speckley domains do not support transport problems");
753 }
754
755 void SpeckleyDomain::setNewX(const escript::Data& arg)
756 {
757 throw SpeckleyException("setNewX(): operation not supported");
758 }
759
760 //protected
761 void SpeckleyDomain::copyData(escript::Data& out, const escript::Data& in) const
762 {
763 const dim_t numComp = in.getDataPointSize();
764 const dim_t numSamples = in.getNumSamples();
765 out.requireWrite();
766 #pragma omp parallel for
767 for (index_t i=0; i<numSamples; i++) {
768 const double* src = in.getSampleDataRO(i);
769 copy(src, src+numComp, out.getSampleDataRW(i));
770 }
771 }
772
773 //protected
774 void SpeckleyDomain::multiplyData(escript::Data& out, const escript::Data& in) const
775 {
776 const dim_t numComp = in.getDataPointSize();
777 const dim_t dpp = out.getNumDataPointsPerSample();
778 const dim_t numSamples = in.getNumSamples();
779 out.requireWrite();
780 #pragma omp parallel for
781 for (index_t i=0; i<numSamples; i++) {
782 const double* src = in.getSampleDataRO(i);
783 double* dest = out.getSampleDataRW(i);
784 for (index_t c=0; c<numComp; c++) {
785 for (index_t q=0; q<dpp; q++)
786 dest[c+q*numComp] = src[c];
787 }
788 }
789 }
790
791 //protected
792 void SpeckleyDomain::updateTagsInUse(int fsType) const
793 {
794 std::vector<int>* tagsInUse=NULL;
795 const std::vector<int>* tags=NULL;
796 switch(fsType) {
797 case Nodes:
798 tags=&m_nodeTags;
799 tagsInUse=&m_nodeTagsInUse;
800 break;
801 case Elements:
802 tags=&m_elementTags;
803 tagsInUse=&m_elementTagsInUse;
804 break;
805 case Points:
806 throw SpeckleyException("updateTagsInUse for Speckley dirac points "
807 "not supported");
808 default:
809 return;
810 }
811
812 // gather global unique tag values from tags into tagsInUse
813 tagsInUse->clear();
814 int lastFoundValue = numeric_limits<int>::min();
815 int minFoundValue, local_minFoundValue;
816 const int numTags = tags->size();
817
818 while (true) {
819 // find smallest value bigger than lastFoundValue
820 minFoundValue = numeric_limits<int>::max();
821 #pragma omp parallel private(local_minFoundValue)
822 {
823 local_minFoundValue = minFoundValue;
824 #pragma omp for schedule(static) nowait
825 for (int i = 0; i < numTags; i++) {
826 const int v = (*tags)[i];
827 if ((v > lastFoundValue) && (v < local_minFoundValue))
828 local_minFoundValue = v;
829 }
830 #pragma omp critical
831 {
832 if (local_minFoundValue < minFoundValue)
833 minFoundValue = local_minFoundValue;
834 }
835 }
836 #ifdef ESYS_MPI
837 local_minFoundValue = minFoundValue;
838 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
839 #endif
840
841 // if we found a new value add it to the tagsInUse vector
842 if (minFoundValue < numeric_limits<int>::max()) {
843 tagsInUse->push_back(minFoundValue);
844 lastFoundValue = minFoundValue;
845 } else
846 break;
847 }
848 }
849
850 //protected
851 void SpeckleyDomain::addToSystemMatrix(escript::AbstractSystemMatrix* mat,
852 const IndexVector& nodes, dim_t numEq,
853 const DoubleVector& array) const
854 {
855 throw SpeckleyException("addToSystemMatrix not yet implemented");
856 }
857
858 //private
859 void SpeckleyDomain::assemblePDE(escript::AbstractSystemMatrix* mat,
860 escript::Data& rhs, const DataMap& coefs,
861 Assembler_ptr assembler) const
862 {
863 if (rhs.isEmpty() && isNotEmpty("X", coefs) && isNotEmpty("Y", coefs))
864 throw SpeckleyException("assemblePDE: right hand side coefficients are "
865 "provided but no right hand side vector given");
866
867 vector<int> fsTypes;
868 assembler->collateFunctionSpaceTypes(fsTypes, coefs);
869
870 if (fsTypes.empty()) {
871 return;
872 }
873
874 int fs=fsTypes[0];
875 if (fs != Elements)
876 throw SpeckleyException("assemblePDE: illegal function space type for coefficients");
877
878 for (vector<int>::const_iterator it=fsTypes.begin()+1; it!=fsTypes.end(); it++) {
879 if (*it != fs) {
880 throw SpeckleyException("assemblePDE: coefficient function spaces don't match");
881 }
882 }
883
884 int numEq, numComp;
885 escript::Data temp(0., rhs.getDataPointShape(), rhs.getFunctionSpace(),
886 rhs.actsExpanded());
887 if (!mat) {
888 if (rhs.isEmpty()) {
889 numEq = numComp = 1;
890 } else {
891 numEq = numComp = rhs.getDataPointSize();
892 }
893 } else {
894 if (!rhs.isEmpty() && rhs.getDataPointSize() != mat->getRowBlockSize())
895 throw SpeckleyException("assemblePDE: matrix row block size and number of components of right hand side don't match");
896 numEq = mat->getRowBlockSize();
897 numComp = mat->getColumnBlockSize();
898 }
899
900 if (numEq != numComp)
901 throw SpeckleyException("assemblePDE: number of equations and number of solutions don't match");
902
903 //TODO: check shape and num samples of coeffs
904
905 if (numEq==1) {
906 assembler->assemblePDESingle(mat, temp, coefs);
907 } else {
908 assembler->assemblePDESystem(mat, temp, coefs);
909 }
910 #ifdef ESYS_MPI
911 balanceNeighbours(temp, false); //summation without averaging
912 #endif
913 rhs += temp; //only now add to rhs because otherwise rhs distorts
914 }
915
916 //private
917 void SpeckleyDomain::assemblePDEBoundary(escript::AbstractSystemMatrix* mat,
918 escript::Data& rhs, const DataMap& coefs, Assembler_ptr assembler) const
919 {
920 if (rhs.isEmpty() && isNotEmpty("y", coefs))
921 throw SpeckleyException("assemblePDEBoundary: y provided but no right hand side vector given");
922
923 int fs=-1;
924 if (isNotEmpty("d", coefs))
925 fs=coefs.find("d")->second.getFunctionSpace().getTypeCode();
926 if (isNotEmpty("y", coefs)) {
927 DataMap::const_iterator iy = coefs.find("y");
928 if (fs == -1)
929 fs = iy->second.getFunctionSpace().getTypeCode();
930 else if (fs != iy->second.getFunctionSpace().getTypeCode())
931 throw SpeckleyException("assemblePDEBoundary: coefficient function spaces don't match");
932 }
933 if (fs==-1) {
934 return;
935 }
936
937 int numEq = 1, numComp = 1;
938 if (!mat) {
939 if (!rhs.isEmpty()) {
940 numEq = numComp = rhs.getDataPointSize();
941 }
942 } else {
943 if (!rhs.isEmpty() && rhs.getDataPointSize() != mat->getRowBlockSize())
944 throw SpeckleyException("assemblePDEBoundary: matrix row block size and number of components of right hand side don't match");
945 numEq = mat->getRowBlockSize();
946 numComp = mat->getColumnBlockSize();
947 }
948
949 if (numEq != numComp)
950 throw SpeckleyException("assemblePDEBoundary: number of equations and number of solutions don't match");
951
952 //TODO: check shape and num samples of coeffs
953
954 if (numEq==1) {
955 assembler->assemblePDEBoundarySingle(mat, rhs, coefs);
956 } else {
957 assembler->assemblePDEBoundarySystem(mat, rhs, coefs);
958 }
959 }
960
961 void SpeckleyDomain::assemblePDEDirac(escript::AbstractSystemMatrix* mat,
962 escript::Data& rhs, const DataMap& coefs,
963 Assembler_ptr assembler) const
964 {
965 bool yNotEmpty = isNotEmpty("y_dirac", coefs);
966 bool dNotEmpty = isNotEmpty("d_dirac", coefs);
967 if (!(yNotEmpty || dNotEmpty)) {
968 return;
969 }
970 escript::Data d = unpackData("d_dirac", coefs);
971 escript::Data y = unpackData("y_dirac", coefs);
972 int nEq = 1, nComp = 1;
973 if (!mat) {
974 if (!rhs.isEmpty()) {
975 nEq = nComp = rhs.getDataPointSize();
976 }
977 } else {
978 if (!rhs.isEmpty() && rhs.getDataPointSize()!=mat->getRowBlockSize())
979 throw SpeckleyException("assemblePDEDirac: matrix row block size "
980 "and number of components of right hand side don't match");
981 nEq = mat->getRowBlockSize();
982 nComp = mat->getColumnBlockSize();
983 }
984
985 for (int i = 0; i < m_diracPoints.size(); i++) { //only for this rank
986 const IndexVector rowIndex(1, m_diracPoints[i].node);
987 if (yNotEmpty) {
988 const double *EM_F = y.getSampleDataRO(i);
989 double *F_p = rhs.getSampleDataRW(0);
990 for (index_t eq = 0; eq < nEq; eq++) {
991 F_p[INDEX2(eq, rowIndex[0], nEq)] += EM_F[INDEX2(eq,i,nEq)];
992 }
993 }
994 if (dNotEmpty) {
995 throw SpeckleyException("Rectangle::assemblePDEDirac currently doesn't support d");
996 // const double *EM_S = d.getSampleDataRO(i);
997 // std::vector<double> contents(EM_S,
998 // EM_S+nEq*nEq*nComp*rowIndex.size()); //TODO: this will break with MPI
999 // addToSystemMatrix(mat, rowIndex, nEq, rowIndex, nComp, contents);
1000 }
1001 }
1002 }
1003
1004
1005
1006 // Expecting ("gaussian", radius, sigma)
1007 bool SpeckleyDomain::supportsFilter(const bp::tuple& t) const
1008 {
1009 if (len(t) == 0) { // so we can handle unfiltered randoms
1010 return true;
1011 }
1012 return false;
1013 }
1014
1015 void SpeckleyDomain::addPoints(const vector<double>& coords,
1016 const vector<int>& tags)
1017 {
1018 for (int i = 0; i < tags.size(); i++) {
1019 dim_t node = findNode(&coords[i * m_numDim]);
1020 if (node >= 0) {
1021 m_diracPointNodeIDs.push_back(borrowSampleReferenceIDs(Nodes)[node]);
1022 DiracPoint dp;
1023 dp.node = node; //local
1024 dp.tag = tags[i];
1025 m_diracPoints.push_back(dp);
1026 } else if (m_mpiInfo->size == 1) {
1027 throw SpeckleyException("Dirac point unmapped, implementation problem in Speckley::addPoints");
1028 }
1029 }
1030 }
1031
1032 } // end of namespace speckley
1033

  ViewVC Help
Powered by ViewVC 1.1.26