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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4211 - (show annotations)
Mon Feb 18 23:54:46 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/plain
File size: 24349 byte(s)
Implemented interpolation from Reduced[Face]Elements to [Face]Elements and
changed regularization to compute gradient on Function instead of
ReducedFunction. Results differ slightly so this should help with the accuracy.

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

  ViewVC Help
Powered by ViewVC 1.1.26