escript Revision_
finley/src/FinleyDomain.h
Go to the documentation of this file.
1
2/*****************************************************************************
3*
4* Copyright (c) 2003-2020 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-2017 by Centre for Geoscience Computing (GeoComp)
14* Development from 2019 by School of Earth and Environmental Sciences
15**
16*****************************************************************************/
17
18#ifndef __FINLEY_DOMAIN_H__
19#define __FINLEY_DOMAIN_H__
20
21/****************************************************************************
22
23 Finley: Domain
24
25 A mesh is built from nodes and elements which describe the domain, surface,
26 and point sources (the latter are needed to establish links with other
27 codes, in particular particle codes). The nodes are stored in a NodeFile
28 and elements in ElementFiles. Finley domains have four ElementFiles
29 containing the elements, surface, contact and point sources, respectively.
30 Notice that the surface elements do not necessarily cover the entire
31 surface of the domain.
32
33 The element type is fixed by the reference element, see ReferenceElement.h.
34 The numbering of the nodes starts with 0.
35
36 Important: it is assumed that every node appears in at least one element or
37 surface element and that any node used in an element, surface element or as
38 a point is specified in the NodeFile, see also resolveNodeIds.
39
40 In some cases it is useful to refer to a mesh entirely built from
41 order 1 (=linear) elements. The linear version of the mesh can be
42 accessed by referring to the first few nodes of each element
43 (thanks to the way the nodes are ordered). As the numbering of
44 these nodes is not continuous a relabeling vector is introduced
45 in the NodeFile. This feature is not fully implemented yet.
46
47 All nodes and elements are tagged. The tag allows to group nodes and
48 elements. A typical application is to mark surface elements on a
49 certain portion of the domain with the same tag. All these surface
50 elements can then be assigned the same value e.g. for the pressure.
51
52 The spatial dimensionality is determined by the type of elements
53 used and can be queried using getDim(). Notice that the element type
54 also determines the type of surface elements to be used.
55
56*****************************************************************************/
57
58#include "system_dep.h"
59
60#include <finley/Finley.h>
61#include <finley/ElementFile.h>
62#include <finley/NodeFile.h>
63#include <finley/Util.h>
64
65#include <escript/AbstractContinuousDomain.h>
66#include <escript/FunctionSpace.h>
67#include <escript/FunctionSpaceFactory.h>
68
69#ifdef ESYS_HAVE_PASO
70#include <paso/SystemMatrixPattern.h>
71#endif
72#ifdef ESYS_HAVE_TRILINOS
73#include <trilinoswrap/types.h>
74#endif
75
76#include <map>
77#include <string>
78#include <vector>
79
80namespace finley {
81
82typedef std::map<std::string, int> TagMap;
83
85 SMT_PASO = 1<<8,
86 SMT_TRILINOS = 1<<10,
87 SMT_COMPLEX = 1<<16,
88 SMT_UNROLL = 1<<17
89};
90
97{
98public:
104 static escript::Domain_ptr load(const std::string& filename);
105
118 static escript::Domain_ptr read(escript::JMPI mpiInfo,
119 const std::string& fileName,
120 int integrationOrder = -1,
121 int reducedIntegrationOrder = -1,
122 bool optimize = false);
123
138 const std::string& filename,
139 int numDim, int integrationOrder = -1,
140 int reducedIntegrationOrder = -1,
141 bool optimize = false,
142 bool useMacroElements = false);
143
161 static escript::Domain_ptr createRec4(dim_t NE0, dim_t NE1,
162 double L0, double L1,
163 bool periodic0, bool periodic1, int order,
164 int reducedOrder, bool useElementsOnFace,
165 bool optimize, escript::JMPI jmpi);
166
187 static escript::Domain_ptr createRec8(dim_t NE0, dim_t NE1,
188 double l0, double l1,
189 bool periodic0, bool periodic1, int order,
190 int reducedOrder, bool useElementsOnFace,
191 bool useFullElementOrder,
192 bool useMacroElements, bool optimize,
193 escript::JMPI jmpi);
194
215 static escript::Domain_ptr createHex8(dim_t NE0, dim_t NE1, dim_t NE2,
216 double l0, double l1, double l2,
217 bool periodic0, bool periodic1, bool periodic2,
218 int order, int reducedOrder,
219 bool useElementsOnFace,
220 bool optimize, escript::JMPI jmpi);
221
244 static escript::Domain_ptr createHex20(dim_t NE0, dim_t NE1, dim_t NE2,
245 double l0, double l1, double l2,
246 bool periodic0, bool periodic1, bool periodic2,
247 int order, int reducedOrder,
248 bool useElementsOnFace,
249 bool useFullElementOrder,
250 bool useMacroElements, bool optimize,
251 escript::JMPI jmpi);
252
261 FinleyDomain(const std::string& name, int numDim, escript::JMPI jmpi);
262
267 FinleyDomain(const FinleyDomain& in);
268
274
280 void addDiracPoints(const std::vector<double>& points,
281 const std::vector<int>& tags);
282
287 NodeFile* getNodes() const { return m_nodes; }
288
293 void setElements(ElementFile* elements);
294
299 ElementFile* getElements() const { return m_elements; }
300
305 void setFaceElements(ElementFile* elements);
306
311 ElementFile* getFaceElements() const { return m_faceElements; }
312
317 void setContactElements(ElementFile* elements);
318
323 ElementFile* getContactElements() const { return m_contactElements; }
324
329 void setPoints(ElementFile* elements);
330
335 ElementFile* getPoints() const { return m_points; }
336
341 virtual escript::JMPI getMPI() const { return m_mpiInfo; }
342
347 virtual int getMPISize() const { return m_mpiInfo->size; }
348
353 virtual int getMPIRank() const { return m_mpiInfo->rank; }
354
359 virtual void MPIBarrier() const;
360
365 virtual bool onMasterProcessor() const { return getMPIRank() == 0; }
366
367 MPI_Comm getMPIComm() const { return m_mpiInfo->comm; }
368
375 void write(const std::string& fileName) const;
376
381 void Print_Mesh_Info(bool full=false) const;
382
388 void dump(const std::string& fileName) const;
389
396 int getTagFromSampleNo(int functionSpaceType, index_t sampleNo) const;
397
403 const index_t* borrowSampleReferenceIDs(int functionSpaceType) const;
404
410 virtual bool isValidFunctionSpaceType(int functionSpaceType) const;
411
416 virtual std::string getDescription() const;
417
422 virtual std::string functionSpaceTypeAsString(int functionSpaceType) const;
423
428 void setFunctionSpaceTypeNames();
429
434 virtual int getContinuousFunctionCode() const;
435
440 virtual int getReducedContinuousFunctionCode() const;
441
446 virtual int getFunctionCode() const;
447
452 virtual int getReducedFunctionCode() const;
453
458 virtual int getFunctionOnBoundaryCode() const;
459
464 virtual int getReducedFunctionOnBoundaryCode() const;
465
470 virtual int getFunctionOnContactZeroCode() const;
471
476 virtual int getReducedFunctionOnContactZeroCode() const;
477
482 virtual int getFunctionOnContactOneCode() const;
483
488 virtual int getReducedFunctionOnContactOneCode() const;
489
494 virtual int getSolutionCode() const;
495
500 virtual int getReducedSolutionCode() const;
501
506 virtual int getDiracDeltaFunctionsCode() const;
507
511 typedef std::map<int, std::string> FunctionSpaceNamesMapType;
512
516 virtual int getDim() const { return m_nodes->numDim; }
517
524 virtual StatusType getStatus() const;
525
530 virtual dim_t getNumDataPointsGlobal() const;
531
537 virtual std::pair<int,dim_t> getDataShape(int functionSpaceCode) const;
538
544 virtual void setToX(escript::Data& arg) const;
545
552 virtual void setTagMap(const std::string& name, int tag);
553
559 virtual int getTag(const std::string& name) const;
560
566 virtual bool isValidTagName(const std::string& name) const;
567
572 virtual std::string showTagNames() const;
573
578 virtual void setNewX(const escript::Data& arg);
579
584 virtual void interpolateOnDomain(escript::Data& target,
585 const escript::Data& source) const;
586
587 virtual bool probeInterpolationOnDomain(int functionSpaceType_source,
588 int functionSpaceType_target) const;
589
590 virtual signed char preferredInterpolationOnDomain(int functionSpaceType_source, int functionSpaceType_target) const;
591
596 bool commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const;
597
602 virtual void interpolateAcross(escript::Data& target, const escript::Data& source) const;
603
607 virtual bool probeInterpolationAcross(int functionSpaceType_source,
608 const escript::AbstractDomain& targetDomain,
609 int functionSpaceType_target) const;
610
616 virtual void setToNormal(escript::Data& out) const;
617
623 virtual void setToSize(escript::Data& out) const;
624
630 virtual void setToGradient(escript::Data& grad, const escript::Data& arg) const;
631
637 virtual void setToIntegrals(std::vector<escript::DataTypes::real_t>& integrals,
638 const escript::Data& arg) const;
639 virtual void setToIntegrals(std::vector<escript::DataTypes::cplx_t>& integrals,
640 const escript::Data& arg) const;
641
650 virtual int getSystemMatrixTypeId(const boost::python::object& options) const;
651
661 virtual int getTransportTypeId(int solver, int preconditioner, int package,
662 bool symmetry) const;
663
669 virtual bool isCellOriented(int functionSpaceCode) const;
670
671 virtual bool ownSample(int fsCode, index_t id) const;
672
677 virtual void addPDEToSystem(
679 const escript::Data& A, const escript::Data& B,
680 const escript::Data& C, const escript::Data& D,
681 const escript::Data& X, const escript::Data& Y,
682 const escript::Data& d, const escript::Data& y,
683 const escript::Data& d_contact,
684 const escript::Data& y_contact,
685 const escript::Data& d_dirac,
686 const escript::Data& y_dirac) const;
687
692 virtual void addPDEToLumpedSystem(escript::Data& mat,
693 const escript::Data& D,
694 const escript::Data& d,
695 const escript::Data& d_dirac,
696 bool useHRZ) const;
697
702 virtual void addPDEToRHS(escript::Data& rhs, const escript::Data& X,
703 const escript::Data& Y, const escript::Data& y,
704 const escript::Data& y_contact,
705 const escript::Data& y_dirac) const;
706
711 virtual void addPDEToTransportProblem(
713 escript::Data& source, const escript::Data& M,
714 const escript::Data& A, const escript::Data& B,
715 const escript::Data& C, const escript::Data& D,
716 const escript::Data& X, const escript::Data& Y,
717 const escript::Data& d, const escript::Data& y,
718 const escript::Data& d_contact,
719 const escript::Data& y_contact,
720 const escript::Data& d_dirac,
721 const escript::Data& y_dirac) const;
722
727 escript::ASM_ptr newSystemMatrix(
728 int row_blocksize,
729 const escript::FunctionSpace& row_functionspace,
730 int column_blocksize,
731 const escript::FunctionSpace& column_functionspace,
732 int type) const;
733
738 escript::ATP_ptr newTransportProblem(int blocksize,
739 const escript::FunctionSpace& functionspace,
740 int type) const;
741
745 virtual escript::Data getX() const;
746
747#ifdef ESYS_HAVE_BOOST_NUMPY
751 virtual boost::python::numpy::ndarray getNumpyX() const;
752
756 virtual boost::python::numpy::ndarray getConnectivityInfo() const;
757#endif
758
762 virtual int getVTKElementType() const;
763
768 virtual escript::Data getNormal() const;
769
773 virtual escript::Data getSize() const;
774
778 virtual bool operator==(const escript::AbstractDomain& other) const;
779 virtual bool operator!=(const escript::AbstractDomain& other) const;
780
785 virtual void setTags(int functionSpaceType, int newTag,
786 const escript::Data& mask) const;
787
793 virtual int getNumberOfTagsInUse(int functionSpaceCode) const;
794
795 virtual const int* borrowListOfTagsInUse(int functionSpaceCode) const;
796
801 virtual bool canTag(int functionSpaceCode) const;
802
806 virtual int getApproximationOrder(int functionSpaceCode) const;
807
808 virtual bool supportsContactElements() const { return true; }
809
810 virtual escript::Data randomFill(const escript::DataTypes::ShapeType& shape,
811 const escript::FunctionSpace& what, long seed,
812 const boost::python::tuple& filter) const;
813
818 const TagMap& getTagMap() const { return m_tagMap; }
819
820 void createMappings(const IndexVector& dofDistribution,
821 const IndexVector& nodeDistribution);
822
823#ifdef ESYS_HAVE_PASO
825 paso::SystemMatrixPattern_ptr getPasoPattern(bool reducedRowOrder,
826 bool reducedColOrder) const;
827#endif
828
829#ifdef ESYS_HAVE_TRILINOS
831 esys_trilinos::const_TrilinosGraph_ptr getTrilinosGraph(bool reducedOrder) const;
832#endif
833
834 void glueFaces(double safetyFactor, double tolerance, bool optimize);
835
836 void joinFaces(double safetyFactor, double tolerance, bool optimize);
837
840 static FinleyDomain* merge(const std::vector<const FinleyDomain*>& meshes);
841
842private:
843 void prepare(bool optimize);
844
845 void setOrders();
846
854 void resolveNodeIds();
855
858 void relabelElementNodes(const IndexVector& newNode, index_t offset);
859
860 template<typename Scalar>
861 void setToIntegralsWorker(std::vector<Scalar>& integrals,
862 const escript::Data& arg) const;
863
864#ifdef ESYS_HAVE_PASO
865 paso::SystemMatrixPattern_ptr makePasoPattern(bool reducedRowOrder,
866 bool reducedColOrder) const;
867#endif
868#ifdef ESYS_HAVE_TRILINOS
869 esys_trilinos::GraphType* createTrilinosGraph(bool reducedOrder) const;
870#endif
871 void createColoring(const IndexVector& dofMap);
872 void distributeByRankOfDOF(const IndexVector& distribution);
873 void markNodes(std::vector<short>& mask, index_t offset, bool useLinear) const;
874 void optimizeDOFDistribution(IndexVector& distribution);
875 void optimizeDOFLabeling(const IndexVector& distribution);
876 void optimizeElementOrdering();
877 void findMatchingFaces(double safetyFactor, double tolerance, int* numPairs,
878 int* elem0, int* elem1, int* matchingNodes) const;
879 void updateTagList();
880 void printElementInfo(const ElementFile* e, const std::string& title,
881 const std::string& defaultType, bool full) const;
882
883 void writeElementInfo(std::ostream& stream, const ElementFile* e,
884 const std::string& defaultType) const;
885
889 std::string m_name;
906#ifdef ESYS_HAVE_PASO
907 // pointer to the sparse matrix patterns
908 mutable paso::SystemMatrixPattern_ptr FullFullPattern;
909 mutable paso::SystemMatrixPattern_ptr FullReducedPattern;
910 mutable paso::SystemMatrixPattern_ptr ReducedFullPattern;
911 mutable paso::SystemMatrixPattern_ptr ReducedReducedPattern;
912#endif
913#ifdef ESYS_HAVE_TRILINOS
914 mutable esys_trilinos::TrilinosGraph_ptr m_fullGraph;
915 mutable esys_trilinos::TrilinosGraph_ptr m_reducedGraph;
916#endif
917
919};
920
921} // end of namespace
922
923#endif // __FINLEY_DOMAIN_H__
924
int MPI_Comm
Definition: EsysMPI.h:44
int getTag(unsigned char sourcex, unsigned char sourcey, unsigned char sourcez, unsigned char targetx, unsigned char targety, unsigned char targetz)
Definition: blocktools.cpp:413
AbstractContinuousDomain, base class for continuous domains.
Definition: AbstractContinuousDomain.h:47
Base class for all escript domains.
Definition: AbstractDomain.h:51
Base class for escript system matrices.
Definition: AbstractSystemMatrix.h:44
Give a short description of what AbstractTransportProblem does.
Definition: AbstractTransportProblem.h:45
Data represents a collection of datapoints.
Definition: Data.h:64
Definition: FunctionSpace.h:36
Definition: finley/src/ElementFile.h:63
FinleyDomain implements the AbstractContinuousDomain interface for the Finley library.
Definition: finley/src/FinleyDomain.h:97
NodeFile * getNodes() const
returns a pointer to this domain's node file
Definition: finley/src/FinleyDomain.h:287
static FunctionSpaceNamesMapType m_functionSpaceTypeNames
Definition: finley/src/FinleyDomain.h:918
int integrationOrder
Definition: finley/src/FinleyDomain.h:892
std::string m_name
domain description
Definition: finley/src/FinleyDomain.h:889
ElementFile * m_faceElements
the table of face elements
Definition: finley/src/FinleyDomain.h:899
NodeFile * m_nodes
the table of the nodes
Definition: finley/src/FinleyDomain.h:895
int reducedApproximationOrder
Definition: finley/src/FinleyDomain.h:891
virtual int getMPIRank() const
returns the number MPI rank of this processor
Definition: finley/src/FinleyDomain.h:353
virtual bool supportsContactElements() const
Definition: finley/src/FinleyDomain.h:808
int approximationOrder
Definition: finley/src/FinleyDomain.h:890
virtual int getMPISize() const
returns the number of processors used for this domain
Definition: finley/src/FinleyDomain.h:347
int reducedIntegrationOrder
Definition: finley/src/FinleyDomain.h:893
ElementFile * getFaceElements() const
returns a pointer to this domain's face element file
Definition: finley/src/FinleyDomain.h:311
virtual int getDim() const
returns the dimensionality of this domain
Definition: finley/src/FinleyDomain.h:516
virtual void setToIntegrals(std::vector< escript::DataTypes::cplx_t > &integrals, const escript::Data &arg) const
ElementFile * m_elements
the table of the elements
Definition: finley/src/FinleyDomain.h:897
ElementFile * getElements() const
returns a pointer to this domain's element file
Definition: finley/src/FinleyDomain.h:299
virtual bool onMasterProcessor() const
returns true if on MPI processor 0, else false
Definition: finley/src/FinleyDomain.h:365
TagMap m_tagMap
the tag map mapping names to tag keys
Definition: finley/src/FinleyDomain.h:905
const TagMap & getTagMap() const
returns a reference to the tag name->value map
Definition: finley/src/FinleyDomain.h:818
ElementFile * m_points
the table of points (treated as elements of dimension 0)
Definition: finley/src/FinleyDomain.h:903
ElementFile * m_contactElements
the table of contact elements
Definition: finley/src/FinleyDomain.h:901
MPI_Comm getMPIComm() const
get the communicator for this domain. Returns an integer on non-MPI builds Routine must be implemente...
Definition: finley/src/FinleyDomain.h:367
std::map< int, std::string > FunctionSpaceNamesMapType
Definition: finley/src/FinleyDomain.h:511
virtual void setToIntegrals(std::vector< escript::DataTypes::real_t > &integrals, const escript::Data &arg) const
copies the integrals of the function defined by arg into integrals. arg has to be defined on this.
virtual escript::JMPI getMPI() const
returns a reference to the MPI information wrapper for this domain
Definition: finley/src/FinleyDomain.h:341
ElementFile * getContactElements() const
returns a pointer to this domain's contact element file
Definition: finley/src/FinleyDomain.h:323
escript::JMPI m_mpiInfo
MPI information.
Definition: finley/src/FinleyDomain.h:887
ElementFile * getPoints() const
returns a pointer to this domain's point (nodal) element file
Definition: finley/src/FinleyDomain.h:335
Definition: finley/src/NodeFile.h:42
#define FINLEY_DLL_API
Definition: finley/src/system_dep.h:29
Domain_ptr readGmsh(const string &fileName, int numDim, int, int, bool optimize)
reads a gmsh mesh file
Definition: dudley/src/DomainFactory.cpp:681
std::vector< index_t > IndexVector
Definition: DataTypes.h:64
std::vector< int > ShapeType
The shape of a single datapoint.
Definition: DataTypes.h:44
index_t dim_t
Definition: DataTypes.h:66
int index_t
type for array/matrix indices used both globally and on each rank
Definition: DataTypes.h:61
Data load(const std::string fileName, const AbstractDomain &domain)
reads Data on domain from file in netCDF format
Definition: DataFactory.cpp:708
boost::shared_ptr< AbstractTransportProblem > ATP_ptr
Definition: AbstractTransportProblem.h:161
boost::shared_ptr< AbstractSystemMatrix > ASM_ptr
Definition: AbstractSystemMatrix.h:35
boost::shared_ptr< AbstractDomain > Domain_ptr
Definition: AbstractDomain.h:43
boost::shared_ptr< JMPI_ > JMPI
Definition: EsysMPI.h:76
A suite of factory methods for creating various finley domains.
Definition: finley/src/Assemble.h:32
SystemMatrixType
Definition: finley/src/FinleyDomain.h:84
@ SMT_TRILINOS
Definition: finley/src/FinleyDomain.h:86
@ SMT_PASO
Definition: finley/src/FinleyDomain.h:85
@ SMT_COMPLEX
Definition: finley/src/FinleyDomain.h:87
@ SMT_UNROLL
Definition: finley/src/FinleyDomain.h:88
std::map< std::string, int > TagMap
Definition: finley/src/FinleyDomain.h:82
Domain_ptr glueFaces(const bp::list &meshList, double safetyFactor, double tolerance, bool optimize)
Definition: finley/src/DomainFactory.cpp:1288
Domain_ptr joinFaces(const bp::list &meshList, double safetyFactor, double tolerance, bool optimize)
Definition: finley/src/DomainFactory.cpp:1300
double l2(dim_t n, const double *x, escript::JMPI mpiinfo)
returns the global L2 norm of x
Definition: PasoUtil.cpp:501
boost::shared_ptr< SystemMatrixPattern > SystemMatrixPattern_ptr
Definition: SystemMatrixPattern.h:41
static dim_t M
Definition: SparseMatrix_saveHB.cpp:37
bool probeInterpolationAcross(int fsType_source, const escript::AbstractDomain &domain, int fsType_target, int dim)
Definition: CrossDomainCoupler.cpp:32