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 
80 namespace finley {
81 
82 typedef 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 {
98 public:
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 
273  ~FinleyDomain();
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 
842 private:
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
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
ElementFile * getContactElements() const
returns a pointer to this domain's contact element file
Definition: finley/src/FinleyDomain.h:323
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
ElementFile * getPoints() const
returns a pointer to this domain's point (nodal) element file
Definition: finley/src/FinleyDomain.h:335
int reducedIntegrationOrder
Definition: finley/src/FinleyDomain.h:893
ElementFile * getElements() const
returns a pointer to this domain's element file
Definition: finley/src/FinleyDomain.h:299
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
virtual bool onMasterProcessor() const
returns true if on MPI processor 0, else false
Definition: finley/src/FinleyDomain.h:365
ElementFile * getFaceElements() const
returns a pointer to this domain's face element file
Definition: finley/src/FinleyDomain.h:311
TagMap m_tagMap
the tag map mapping names to tag keys
Definition: finley/src/FinleyDomain.h:905
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
NodeFile * getNodes() const
returns a pointer to this domain's node file
Definition: finley/src/FinleyDomain.h:287
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
const TagMap & getTagMap() const
returns a reference to the tag name->value map
Definition: finley/src/FinleyDomain.h:818
escript::JMPI m_mpiInfo
MPI information.
Definition: finley/src/FinleyDomain.h:887
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:33
boost::shared_ptr< AbstractDomain > Domain_ptr
Definition: AbstractDomain.h:41
boost::shared_ptr< JMPI_ > JMPI
Definition: EsysMPI.h:74
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:40
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