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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4622 - (show annotations)
Fri Jan 17 04:55:41 2014 UTC (5 years, 3 months ago) by sshaw
File MIME type: text/plain
File size: 24912 byte(s)
Added dirac support to ripley, added interface for custom assemblers for ripleydomains (also added the custom assembler for 2D VTI waves), changed synthetic_VTI example to use the new, faster custom assembler

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

  ViewVC Help
Powered by ViewVC 1.1.26