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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6141 - (show annotations)
Wed Apr 6 03:51:30 2016 UTC (2 years, 8 months ago) by caltinay
File MIME type: text/plain
File size: 24877 byte(s)
more namespacing of defines.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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 #ifndef __Speckley_DOMAIN_H__
18 #define __Speckley_DOMAIN_H__
19
20 #include <speckley/Speckley.h>
21 #include <speckley/SpeckleyException.h>
22 #include <speckley/AbstractAssembler.h>
23 #include <speckley/domainhelpers.h>
24
25 #include <escript/AbstractContinuousDomain.h>
26 #include <escript/Data.h>
27 #include <escript/FunctionSpace.h>
28 #include <escript/SubWorld.h>
29
30 #include <boost/python/tuple.hpp>
31 #include <boost/python/list.hpp>
32
33 namespace speckley {
34
35 enum assembler_t {
36 DEFAULT_ASSEMBLER
37 };
38
39 /* There is no particular significance to this type,
40 It is here as a typedef because a bug in clang++ prevents
41 that compiler from recognising it as a valid part of
42 a constant expression.
43 */
44 typedef std::map<std::string, int> simap_t;
45
46
47 /**
48 \brief
49 Structure that wraps parameters for the grid reading routines.
50 */
51 struct ReaderParameters
52 {
53 /// the (global) offset into the data object to start writing into
54 std::vector<dim_t> first;
55 /// the number of values to read from file
56 std::vector<dim_t> numValues;
57 /// how often to write each value from the file into the data object
58 /// (e.g. to supersample)
59 std::vector<int> multiplier;
60 /// if non-zero, values are written from last index to first index
61 std::vector<int> reverse;
62 /// byte order in the file (used by binary reader only)
63 int byteOrder;
64 /// data type in the file (used by binary reader only)
65 int dataType;
66 };
67
68 /**
69 \brief
70 A struct to contain a dirac point's information.
71 */
72 struct DiracPoint
73 {
74 dim_t node;
75 int tag;
76 };
77
78 /**
79 \brief
80 SpeckleyDomain extends the AbstractContinuousDomain interface
81 for the Speckley library and is the base class for Rectangle and Brick.
82 */
83
84 class Speckley_DLL_API SpeckleyDomain : public escript::AbstractContinuousDomain
85 {
86 public:
87 /**
88 \brief
89 Constructor with number of dimensions. Allocates MPI info structure.
90 */
91 SpeckleyDomain(dim_t dim, int order, escript::SubWorld_ptr p=escript::SubWorld_ptr());
92
93 /**
94 \brief
95 Destructor
96 */
97 ~SpeckleyDomain();
98
99 /**
100 \brief
101 returns a reference to the MPI information wrapper for this domain
102 */
103 virtual escript::JMPI getMPI() const { return m_mpiInfo; }
104
105 /**
106 \brief
107 returns the number of processors used for this domain
108 */
109 virtual int getMPISize() const { return m_mpiInfo->size; }
110
111 /**
112 \brief
113 returns the MPI rank of this processor
114 */
115 virtual int getMPIRank() const { return m_mpiInfo->rank; }
116
117 /**
118 \brief
119 if compiled for MPI then executes an MPI_Barrier, else does nothing
120 */
121 virtual void MPIBarrier() const {
122 #ifdef ESYS_MPI
123 MPI_Barrier(m_mpiInfo->comm);
124 #endif
125 }
126
127 /**
128 \brief
129 returns true if on MPI processor 0, else false
130 */
131 virtual bool onMasterProcessor() const { return getMPIRank()==0; }
132
133 /**
134 \brief
135 returns the MPI communicator
136 */
137 MPI_Comm getMPIComm() const { return m_mpiInfo->comm; }
138
139 /**
140 \brief
141 returns true if the argument is a valid function space type for this
142 domain
143 */
144 virtual bool isValidFunctionSpaceType(int fsType) const;
145
146 /**
147 \brief
148 returns a description for the given function space type code
149 */
150 virtual std::string functionSpaceTypeAsString(int fsType) const;
151
152 /**
153 \brief
154 returns the number of spatial dimensions of the domain
155 */
156 virtual int getDim() const { return m_numDim; }
157
158 /**
159 \brief equality operator
160 */
161 virtual bool operator==(const escript::AbstractDomain& other) const;
162
163 /**
164 \brief inequality operator
165 */
166 virtual bool operator!=(const escript::AbstractDomain& other) const {
167 return !(operator==(other));
168 }
169
170 /**
171 \brief
172 returns the number of data points per sample, and the number of samples
173 as a pair.
174 \param fsType The function space type
175 */
176 virtual std::pair<int,dim_t> getDataShape(int fsType) const;
177
178 /**
179 \brief
180 returns the tag key for the given sample number
181 \param fsType The function space type
182 \param sampleNo The sample number
183 */
184 int getTagFromSampleNo(int fsType, dim_t sampleNo) const;
185
186 /**
187 \brief
188 sets a map from a clear tag name to a tag key
189 \param name tag name
190 \param tag tag key
191 */
192 virtual void setTagMap(const std::string& name, int tag) {
193 m_tagMap[name] = tag;
194 }
195
196 /**
197 \brief
198 returns the tag key for tag name
199 \param name tag name
200 */
201 virtual int getTag(const std::string& name) const {
202 if (m_tagMap.find(name) != m_tagMap.end()) {
203 return m_tagMap.find(name)->second;
204 } else {
205 throw SpeckleyException("getTag: invalid tag name");
206 }
207 }
208
209 /**
210 \brief
211 returns true if name is a defined tag name
212 \param name tag name to be checked
213 */
214 virtual bool isValidTagName(const std::string& name) const {
215 return (m_tagMap.find(name)!=m_tagMap.end());
216 }
217
218 /**
219 \brief
220 returns all tag names in a single string separated by commas
221 */
222 virtual std::string showTagNames() const;
223
224 /**
225 \brief
226 assigns new location to the domain.
227 \note This is not supported in Speckley
228 */
229 virtual void setNewX(const escript::Data& arg);
230
231 /**
232 \brief
233 interpolates data given on source onto target where source and target
234 have to be given on the same domain
235 */
236 virtual void interpolateOnDomain(escript::Data& target,
237 const escript::Data& source) const;
238
239 /**
240 \brief
241 returns true if data on fsType_source can be interpolated onto
242 fsType_target, false otherwise
243 */
244 virtual bool probeInterpolationOnDomain(int fsType_source,
245 int fsType_target) const;
246
247 /**
248 \brief Preferred direction of interpolation.
249 If you really need to test for a particular direction, then use
250 probeInterpolation.
251 \return 0 for not possible, 1 for possible and preferred, -1 other
252 direction preferred (does not mean this direction is possible)
253 */
254 virtual signed char preferredInterpolationOnDomain(int fsType_source,
255 int fsType_target) const;
256
257 /**
258 \brief
259 given a vector of FunctionSpace type codes, passes back a code which all
260 can be interpolated to
261 \return true if result is valid, false if not
262 */
263 bool
264 commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const;
265
266 /**
267 \brief
268 interpolates data given on source onto target where source and target
269 are given on different domains
270 */
271 virtual void interpolateAcross(escript::Data& target,
272 const escript::Data& source) const = 0;
273
274 /**
275 \brief
276 determines whether interpolation from source to target is possible
277 */
278 virtual bool probeInterpolationAcross(int, const escript::AbstractDomain&,
279 int) const = 0;
280
281 /**
282 \brief
283 returns locations in the SEM nodes
284 */
285 virtual escript::Data getX() const;
286
287 /**
288 \brief
289 returns boundary normals at the quadrature point on the face elements
290 */
291 virtual escript::Data getNormal() const;
292
293 /**
294 \brief returns the element size
295 */
296 virtual escript::Data getSize() const;
297
298 /**
299 \brief
300 copies the location of data points into arg. The domain of arg has to
301 match this domain.
302 */
303 virtual void setToX(escript::Data& arg) const;
304
305 /**
306 \brief
307 copies the gradient of 'in' into 'out'. The actual function space to be
308 considered for the gradient is defined by 'in'. Both arguments have to
309 be defined on this domain.
310 */
311 virtual void setToGradient(escript::Data& out,
312 const escript::Data& in) const;
313
314 /**
315 \brief
316 assigns new tag newTag to all samples of given function space with a
317 positive value of mask for any of its sample points
318 */
319 virtual void setTags(int fsType, int newTag, const escript::Data& mask) const;
320
321 /**
322 \brief
323 returns true if data on this domain and given function space type has
324 to be considered as cell centered data
325 */
326 virtual bool isCellOriented(int fsType) const;
327
328 /**
329 \brief
330 returns a status indicator of the domain. The status identifier should
331 be unique over the lifetime of the object but may be updated if changes
332 to the domain happen, e.g. modifications to its geometry.
333 */
334 virtual StatusType getStatus() const { return m_status; }
335
336 /**
337 \brief
338 returns the number of tags in use for a function space type
339 */
340 virtual int getNumberOfTagsInUse(int fsType) const;
341
342 /**
343 \brief
344 returns a pointer to the list of tags in use for a function space type
345 */
346 virtual const int* borrowListOfTagsInUse(int fsType) const;
347
348 /**
349 \brief
350 checks if this domain allows tags for the specified function space type
351 */
352 virtual bool canTag(int fsType) const;
353
354 /**
355 \brief
356 returns the approximation order used for a function space
357 */
358 virtual int getApproximationOrder(int fsType) const { return 1; }
359
360 /**
361 \brief
362 returns true if this domain supports contact elements, false otherwise
363 */
364 virtual bool supportsContactElements() const { return false; }
365
366 /**
367 \brief
368 returns a continuous FunctionSpace code
369 */
370 virtual int getContinuousFunctionCode() const { return Nodes; }
371
372 /**
373 \brief
374 returns a continuous on reduced order nodes FunctionSpace code
375 */
376 virtual int getReducedContinuousFunctionCode() const {
377 throw SpeckleyException("Speckley does not support reduced functionspaces");
378 }
379
380 /**
381 \brief
382 returns a function FunctionSpace code
383 */
384 virtual int getFunctionCode() const { return Elements; }
385
386 /**
387 \brief
388 returns a function with reduced integration order FunctionSpace code
389 */
390 virtual int getReducedFunctionCode() const {
391 return ReducedElements;
392 }
393
394 /**
395 \brief
396 returns a function on boundary FunctionSpace code
397 */
398 virtual int getFunctionOnBoundaryCode() const {
399 throw SpeckleyException("Speckley does not support face elements");
400 }
401
402 /**
403 \brief
404 returns a function on boundary with reduced integration order
405 FunctionSpace code
406 */
407 virtual int getReducedFunctionOnBoundaryCode() const {
408 throw SpeckleyException("Speckley does not support face elements");
409 }
410
411 /**
412 \brief
413 return a FunctionOnContactZero code
414 */
415 virtual int getFunctionOnContactZeroCode() const {
416 throw SpeckleyException("Speckley does not support contact elements");
417 }
418
419 /**
420 \brief
421 returns a FunctionOnContactZero code with reduced integration order
422 */
423 virtual int getReducedFunctionOnContactZeroCode() const {
424 throw SpeckleyException("Speckley does not support contact elements");
425 }
426
427 /**
428 \brief
429 returns a FunctionOnContactOne code
430 */
431 virtual int getFunctionOnContactOneCode() const {
432 throw SpeckleyException("Speckley does not support contact elements");
433 }
434
435 /**
436 \brief
437 returns a FunctionOnContactOne code with reduced integration order
438 */
439 virtual int getReducedFunctionOnContactOneCode() const {
440 throw SpeckleyException("Speckley does not support contact elements");
441 }
442
443 /**
444 \brief
445 returns a Solution FunctionSpace code
446 */
447 virtual int getSolutionCode() const {
448 return DegreesOfFreedom;
449 }
450
451 /**
452 \brief
453 returns a ReducedSolution FunctionSpace code
454 */
455 virtual int getReducedSolutionCode() const {
456 throw SpeckleyException("Speckley does not support reduced function spaces");
457 }
458
459 /**
460 \brief
461 returns a DiracDeltaFunctions FunctionSpace code
462 */
463 virtual int getDiracDeltaFunctionsCode() const { return Points; }
464
465 /**
466 \brief
467 returns the identifier of the matrix type to be used for the global
468 stiffness matrix when a particular solver, package, preconditioner,
469 and symmetric matrix is used
470 \param options a python object containing the solver, package,
471 preconditioner and symmetry
472 */
473 virtual int getSystemMatrixTypeId(const boost::python::object& options) const;
474
475 /**
476 \brief
477 returns the identifier of the transport problem type to be used when a
478 particular solver, preconditioner, package and symmetric matrix is used
479 \param solver
480 \param preconditioner
481 \param package
482 \param symmetry
483 */
484 virtual int getTransportTypeId(int solver, int preconditioner, int package,
485 bool symmetry) const;
486
487 /**
488 \brief
489 copies the integrals of the function defined by arg into integrals.
490 arg has to be defined on this domain.
491 */
492 virtual void setToIntegrals(DoubleVector& integrals, const escript::Data& arg) const;
493
494
495 /**
496 \brief
497 adds a PDE onto the stiffness matrix mat and rhs, used for custom
498 solvers with varying arguments counts and so on
499 */
500 virtual void addToSystem(escript::AbstractSystemMatrix& mat,
501 escript::Data& rhs, const DataMap& data,
502 Assembler_ptr assembler) const;
503
504 /**
505 \brief
506 a wrapper for addToSystem that allows calling from Python
507 */
508 virtual void addToSystemFromPython(escript::AbstractSystemMatrix& mat,
509 escript::Data& rhs, const boost::python::list& data,
510 Assembler_ptr assembler) const;
511
512 /**
513 \brief
514 adds a PDE onto rhs, used for custom
515 solvers with varying arguments counts and so on
516 */
517 virtual void addToRHS(escript::Data& rhs, const DataMap& data,
518 Assembler_ptr assembler) const;
519
520 /**
521 \brief
522 a wrapper for addToRHS that allows calling from Python
523 */
524 virtual void addToRHSFromPython(escript::Data& rhs,
525 const boost::python::list& data,
526 Assembler_ptr assembler) const;
527
528 using escript::AbstractContinuousDomain::addPDEToTransportProblem;
529 /**
530 \brief
531 adds a PDE onto a transport problem
532 */
533 virtual void addPDEToTransportProblem(escript::AbstractTransportProblem& tp,
534 escript::Data& source, const DataMap& data,
535 Assembler_ptr assembler) const;
536 /**
537 \brief
538 adds a PDE onto a transport problem
539 */
540 void addPDEToTransportProblemFromPython(
541 escript::AbstractTransportProblem& tp,
542 escript::Data& source, const boost::python::list& data,
543 Assembler_ptr assembler) const;
544
545 /**
546 \brief
547 creates a stiffness matrix and initializes it with zeros
548 */
549 virtual escript::ASM_ptr newSystemMatrix(int row_blocksize,
550 const escript::FunctionSpace& row_functionspace,
551 int column_blocksize,
552 const escript::FunctionSpace& column_functionspace, int type) const;
553
554 /**
555 \brief
556 creates a transport problem
557 */
558 virtual escript::ATP_ptr newTransportProblem(int blocksize,
559 const escript::FunctionSpace& functionspace,
560 int type) const;
561
562 /**
563 \brief
564 writes information about the mesh to standard output
565 \param full whether to print additional data
566 */
567 virtual void Print_Mesh_Info(bool full=false) const;
568
569
570 /************************************************************************/
571
572 /**
573 \brief
574 writes the current mesh to a file with the given name
575 \param filename The name of the file to write to
576 */
577 virtual void write(const std::string& filename) const = 0;
578
579 /**
580 \brief
581 returns a description for this domain
582 */
583 virtual std::string getDescription() const = 0;
584
585 /**
586 \brief
587 dumps the mesh to a file with the given name
588 \param filename The name of the output file
589 */
590 void dump(const std::string& filename) const = 0;
591
592 /**
593 \brief
594 returns the array of reference numbers for a function space type
595 \param fsType The function space type
596 */
597 const index_t* borrowSampleReferenceIDs(int fsType) const = 0;
598
599 /**
600 \brief
601 copies the surface normals at data points into out. The actual function
602 space to be considered is defined by out. out has to be defined on this
603 domain.
604 */
605 virtual void setToNormal(escript::Data& out) const = 0;
606
607 /**
608 \brief
609 copies the size of samples into out. The actual function space to be
610 considered is defined by out. out has to be defined on this domain.
611 */
612 virtual void setToSize(escript::Data& out) const = 0;
613
614 /**
615 \brief
616 reads grid data from a netCDF file into a Data object
617 */
618 virtual void readNcGrid(escript::Data& out, std::string filename,
619 std::string varname, const ReaderParameters& params) const = 0;
620
621 /**
622 \brief
623 reads grid data from a raw binary file into a Data object
624 */
625 virtual void readBinaryGrid(escript::Data& out, std::string filename,
626 const ReaderParameters& params) const = 0;
627
628 /**
629 \brief
630 reads grid data from a compressed raw binary file into a Data object
631 */
632 virtual void readBinaryGridFromZipped(escript::Data& out,
633 std::string filename, const ReaderParameters& params) const = 0;
634
635 /**
636 \brief
637 writes a Data object to a file in raw binary format
638 */
639 virtual void writeBinaryGrid(const escript::Data& in, std::string filename,
640 int byteOrder, int dataType) const = 0;
641
642 /**
643 \brief
644 returns true if this rank owns the sample id on given function space
645 */
646 virtual bool ownSample(int fsType, index_t id) const = 0;
647
648 /**
649 \brief
650 returns the number of data points summed across all MPI processes
651 */
652 virtual dim_t getNumDataPointsGlobal() const = 0;
653
654 /**
655 \brief
656 returns the number of nodes per MPI rank in each dimension
657 */
658 virtual const dim_t* getNumNodesPerDim() const = 0;
659
660 /**
661 \brief
662 returns the number of elements per MPI rank in each dimension
663 */
664 virtual const dim_t* getNumElementsPerDim() const = 0;
665
666 /**
667 \brief
668 returns the number of face elements in the order
669 (left,right,bottom,top,[front,back]) on current MPI rank
670 */
671 virtual const dim_t* getNumFacesPerBoundary() const = 0;
672
673 /**
674 \brief
675 returns the node distribution vector
676 */
677 virtual IndexVector getNodeDistribution() const = 0;
678
679 /**
680 \brief
681 returns the number of spatial subdivisions in each dimension
682 */
683 virtual const int* getNumSubdivisionsPerDim() const = 0;
684
685 /**
686 \brief
687 returns the index'th coordinate value in given dimension for this rank
688 */
689 virtual double getLocalCoordinate(dim_t index, int dim) const = 0;
690
691 /**
692 \brief
693 returns the tuple (origin, spacing, number_of_elements)
694 */
695 virtual boost::python::tuple getGridParameters() const = 0;
696
697 /**
698 \brief
699 returns true if this domain can handle the specified tuple of filter
700 options.
701 */
702 virtual bool supportsFilter(const boost::python::tuple& t) const;
703
704 /**
705 \brief
706 */
707 virtual Assembler_ptr createAssembler(const std::string type,
708 const DataMap& options) const {
709 throw SpeckleyException("Domain does not support custom assemblers");
710 }
711
712 Assembler_ptr createAssemblerFromPython(const std::string type,
713 const boost::python::list& options) const;
714
715 /**
716 \brief
717 returns the lengths of the domain
718 */
719 virtual const double *getLength() const = 0;
720
721 /**
722 \brief
723 returns the order of the domain
724 */
725 inline int getOrder() const { return m_order;}
726
727 protected:
728 int m_numDim;
729 StatusType m_status;
730 escript::JMPI m_mpiInfo;
731 TagMap m_tagMap;
732 mutable std::vector<int> m_nodeTags, m_nodeTagsInUse;
733 mutable std::vector<int> m_elementTags, m_elementTagsInUse;
734 std::vector<DiracPoint> m_diracPoints;
735 IndexVector m_diracPointNodeIDs; //for borrowSampleID
736 assembler_t assembler_type;
737 /// element order (will be m_order + 1 quad points in each axis)
738 int m_order;
739
740 /// copies data in 'in' to 'out' (both must be on same function space)
741 void copyData(escript::Data& out, const escript::Data& in) const;
742
743 // this is const because setTags is const
744 void updateTagsInUse(int fsType) const;
745
746 void addToSystemMatrix(escript::AbstractSystemMatrix* mat,
747 const IndexVector& nodes, dim_t numEq,
748 const DoubleVector& array) const;
749
750 void addPoints(const std::vector<double>& coords,
751 const std::vector<int>& tags);
752
753 /// expands ReducedFunction (in) to Function (out)
754 void multiplyData(escript::Data& out, const escript::Data& in) const;
755
756 /***********************************************************************/
757
758 /// returns the number of nodes per MPI rank
759 virtual dim_t getNumNodes() const = 0;
760
761 /// returns the number of elements per MPI rank
762 virtual dim_t getNumElements() const = 0;
763
764 /// returns the number of degrees of freedom per MPI rank
765 virtual dim_t getNumDOF() const = 0;
766
767 /// populates the data object 'arg' with the node coordinates
768 virtual void assembleCoordinates(escript::Data& arg) const = 0;
769
770 /// computes the gradient of 'in' and puts the result in 'out'
771 virtual void assembleGradient(escript::Data& out,
772 const escript::Data& in) const = 0;
773
774 /// copies the integrals of the function defined by 'arg' into 'integrals'
775 virtual void assembleIntegrate(DoubleVector& integrals,
776 const escript::Data& arg) const = 0;
777
778 /// interpolates data on nodes in 'in' onto elements in 'out'
779 virtual void interpolateNodesOnElements(escript::Data& out,
780 const escript::Data& in,
781 bool reduced) const = 0;
782
783 /// interpolates data on elements in 'in' onto nodes in 'out'
784 virtual void interpolateElementsOnNodes(escript::Data& out,
785 const escript::Data& in) const = 0;
786
787 virtual dim_t getDofOfNode(dim_t node) const = 0;
788
789 /// interpolates from Element -> ReducedElement
790 virtual void reduceElements(escript::Data& out, const escript::Data& in) const = 0;
791
792 #ifdef ESYS_MPI
793 /// sum up overlapping edges of MPI ranks and average if average is true
794 virtual void balanceNeighbours(escript::Data& data, bool average) const = 0;
795 #endif
796
797 private:
798 /// calls the right PDE assembly routines after performing input checks
799 void assemblePDE(escript::AbstractSystemMatrix* mat, escript::Data& rhs,
800 const DataMap& coefs, Assembler_ptr assembler) const;
801
802 /// calls the right PDE boundary assembly routines after performing input
803 /// checks
804 void assemblePDEBoundary(escript::AbstractSystemMatrix* mat,
805 escript::Data& rhs, const DataMap& coefs,
806 Assembler_ptr assembler) const;
807
808 void assemblePDEDirac(escript::AbstractSystemMatrix* mat,
809 escript::Data& rhs, const DataMap& coefs,
810 Assembler_ptr assembler) const;
811
812 /// finds the node that the given point coordinates belong to
813 virtual dim_t findNode(const double *coords) const = 0;
814 };
815
816 } // end of namespace speckley
817
818 #endif // __Speckley_DOMAIN_H__
819

  ViewVC Help
Powered by ViewVC 1.1.26