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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5478 - (show annotations)
Wed Feb 18 01:57:49 2015 UTC (4 years, 6 months ago) by sshaw
File MIME type: text/plain
File size: 13483 byte(s)
introducing shiny new reduced function spaces to speckley
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2015 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 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #ifndef __Speckley_BRICK_H__
18 #define __Speckley_BRICK_H__
19
20 #include <speckley/SpeckleyDomain.h>
21
22 namespace speckley {
23
24 #ifdef USE_RIPLEY
25 class RipleyCoupler; //forward declaration of coupler to avoid circles
26 #endif
27
28 /**
29 \brief
30 Brick is the 3-dimensional implementation of a SpeckleyDomain.
31 */
32 class Speckley_DLL_API Brick: public SpeckleyDomain
33 {
34 friend class DefaultAssembler3D;
35 public:
36
37 /**
38 \brief creates a hexagonal mesh with n0 x n1 x n2 elements over the
39 brick [x0,x1] x [y0,y1] x [z0,z1].
40 \param n0,n1,n2 number of elements in each dimension
41 \param x0,y0,z0,x1,y1,z1 coordinates of corner nodes of the brick
42 \param d0,d1,d2 number of subdivisions in each dimension
43 */
44 Brick(int order, dim_t n0, dim_t n1, dim_t n2, double x0, double y0, double z0, double x1,
45 double y1, double z1, int d0=-1, int d1=-1, int d2=-1,
46 const std::vector<double>& points = std::vector<double>(),
47 const std::vector<int>& tags = std::vector<int>(),
48 const TagMap& tagnamestonums = TagMap(),
49 escript::SubWorld_ptr w=escript::SubWorld_ptr());
50
51 /**
52 \brief
53 Destructor.
54 */
55 ~Brick();
56
57 /**
58 \brief
59 returns a description for this domain
60 */
61 virtual std::string getDescription() const;
62
63 /**
64 \brief equality operator
65 */
66 virtual bool operator==(const escript::AbstractDomain& other) const;
67
68 /**
69 \brief
70 writes the current mesh to a file with the given name
71 \param filename The name of the file to write to
72 */
73 virtual void write(const std::string& filename) const;
74
75 /**
76 \brief
77 dumps the mesh to a file with the given name
78 \param filename The name of the output file
79 */
80 void dump(const std::string& filename) const;
81
82 /**
83 */
84 virtual void readNcGrid(escript::Data& out, std::string filename,
85 std::string varname, const ReaderParameters& params) const;
86
87 /**
88 */
89 virtual void readBinaryGrid(escript::Data& out, std::string filename,
90 const ReaderParameters& params) const;
91
92 #ifdef USE_BOOSTIO
93 virtual void readBinaryGridFromZipped(escript::Data& out, std::string filename,
94 const ReaderParameters& params) const;
95 #endif
96
97 /**
98 */
99 virtual void writeBinaryGrid(const escript::Data& in,
100 std::string filename,
101 int byteOrder, int dataType) const;
102
103 /**
104 \brief
105 returns the array of reference numbers for a function space type
106 \param fsType The function space type
107 */
108 const dim_t* borrowSampleReferenceIDs(int fsType) const;
109
110 /**
111 \brief
112 returns true if this rank owns the sample id.
113 */
114 virtual bool ownSample(int fsType, index_t id) const;
115
116 /**
117 \brief
118 copies the surface normals at data points into out. The actual function
119 space to be considered is defined by out. out has to be defined on this
120 domain.
121 */
122 virtual void setToNormal(escript::Data& out) const;
123
124 /**
125 \brief
126 copies the size of samples into out. The actual function space to be
127 considered is defined by out. out has to be defined on this domain.
128 */
129 virtual void setToSize(escript::Data& out) const;
130
131 /**
132 \brief
133 returns the number of data points summed across all MPI processes
134 */
135 virtual dim_t getNumDataPointsGlobal() const;
136
137 /**
138 \brief
139 writes information about the mesh to standard output
140 \param full whether to print additional data
141 */
142 virtual void Print_Mesh_Info(const bool full=false) const;
143
144 /**
145 \brief
146 returns the number of nodes per MPI rank in each dimension
147 */
148 virtual const dim_t* getNumNodesPerDim() const { return m_NN; }
149
150 /**
151 \brief
152 returns the number of elements per MPI rank in each dimension
153 */
154 virtual const dim_t* getNumElementsPerDim() const { return m_NE; }
155
156 /**
157 \brief
158 returns the number of face elements in the order
159 (left,right,bottom,top,front,back) on current MPI rank
160 */
161 virtual const dim_t* getNumFacesPerBoundary() const { return m_faceCount; }
162
163 /**
164 \brief
165 returns the node distribution vector
166 */
167 virtual IndexVector getNodeDistribution() const { return m_nodeDistribution; }
168
169 /**
170 \brief
171 returns the number of spatial subdivisions in each dimension
172 */
173 virtual const int* getNumSubdivisionsPerDim() const { return m_NX; }
174
175 /**
176 \brief
177 returns the index'th coordinate value in given dimension for this rank
178 */
179 virtual double getLocalCoordinate(index_t index, int dim) const;
180
181 /**
182 \brief
183 returns the tuple (origin, spacing, number_of_elements)
184 */
185 virtual boost::python::tuple getGridParameters() const;
186
187 /**
188 * \brief
189 Returns a Data object filled with random data passed through filter.
190 */
191 virtual escript::Data randomFill(const escript::DataTypes::ShapeType& shape,
192 const escript::FunctionSpace& what,
193 long seed,
194 const boost::python::tuple& filter) const;
195
196 /**
197 \brief
198 interpolates data given on source onto target where source and target
199 are given on different domains
200 */
201 virtual void interpolateAcross(escript::Data& target,
202 const escript::Data& source) const;
203
204 /**
205 \brief
206 determines whether interpolation from source to target is possible
207 */
208 virtual bool probeInterpolationAcross(int, const escript::AbstractDomain&,
209 int) const;
210
211 /**
212 \brief
213 returns the lengths of the domain
214 */
215 const double *getLength() const { return m_length; }
216
217 protected:
218 virtual dim_t getNumNodes() const;
219 virtual dim_t getNumElements() const;
220 virtual dim_t getNumDOF() const;
221 virtual void assembleCoordinates(escript::Data& arg) const;
222 virtual void assembleGradient(escript::Data& out,
223 const escript::Data& in) const;
224 virtual void assembleIntegrate(DoubleVector& integrals,
225 const escript::Data& arg) const;
226 virtual void interpolateNodesOnElements(escript::Data& out,
227 const escript::Data& in,
228 bool reduced) const;
229 virtual void interpolateElementsOnNodes(escript::Data& out,
230 const escript::Data& in) const;
231 virtual dim_t getDofOfNode(dim_t node) const;
232 Assembler_ptr createAssembler(std::string type, const DataMap& constants) const;
233 virtual void reduceElements(escript::Data& out, const escript::Data& in) const;
234 #ifdef ESYS_MPI
235 virtual void balanceNeighbours(escript::Data& data, bool average) const;
236 #endif
237
238 private:
239 void gradient_order2(escript::Data&, const escript::Data&) const;
240 void gradient_order3(escript::Data&, const escript::Data&) const;
241 void gradient_order4(escript::Data&, const escript::Data&) const;
242 void gradient_order5(escript::Data&, const escript::Data&) const;
243 void gradient_order6(escript::Data&, const escript::Data&) const;
244 void gradient_order7(escript::Data&, const escript::Data&) const;
245 void gradient_order8(escript::Data&, const escript::Data&) const;
246 void gradient_order9(escript::Data&, const escript::Data&) const;
247 void gradient_order10(escript::Data&, const escript::Data&) const;
248
249 void reduction_order2(const escript::Data&, escript::Data&) const;
250 void reduction_order3(const escript::Data&, escript::Data&) const;
251 void reduction_order4(const escript::Data&, escript::Data&) const;
252 void reduction_order5(const escript::Data&, escript::Data&) const;
253 void reduction_order6(const escript::Data&, escript::Data&) const;
254 void reduction_order7(const escript::Data&, escript::Data&) const;
255 void reduction_order8(const escript::Data&, escript::Data&) const;
256 void reduction_order9(const escript::Data&, escript::Data&) const;
257 void reduction_order10(const escript::Data&, escript::Data&) const;
258
259 void integral_order2(std::vector<double>&, const escript::Data&) const;
260 void integral_order3(std::vector<double>&, const escript::Data&) const;
261 void integral_order4(std::vector<double>&, const escript::Data&) const;
262 void integral_order5(std::vector<double>&, const escript::Data&) const;
263 void integral_order6(std::vector<double>&, const escript::Data&) const;
264 void integral_order7(std::vector<double>&, const escript::Data&) const;
265 void integral_order8(std::vector<double>&, const escript::Data&) const;
266 void integral_order9(std::vector<double>&, const escript::Data&) const;
267 void integral_order10(std::vector<double>&, const escript::Data&) const;
268 #ifdef ESYS_MPI
269 void setCornerNeighbours();
270 void shareEdges(escript::Data& out, int rx, int ry, int rz) const;
271 void shareFaces(escript::Data& out, int rx, int ry, int rz) const;
272 void shareCorners(escript::Data& out) const;
273 #endif
274 /* \brief
275 interpolates the non-corner point values of an element
276 from the corner values
277 */
278 void interpolateFromCorners(escript::Data& out) const;
279
280 void populateSampleIds();
281 void addToMatrixAndRHS(escript::AbstractSystemMatrix* S, escript::Data& F,
282 const DoubleVector& EM_S, const DoubleVector& EM_F,
283 bool addS, bool addF, index_t firstNode, int nEq=1, int nComp=1) const;
284
285 template<typename ValueType>
286 void readBinaryGridImpl(escript::Data& out, const std::string& filename,
287 const ReaderParameters& params) const;
288 template<typename ValueType>
289 void readBinaryGridZippedImpl(escript::Data& out,
290 const std::string& filename, const ReaderParameters& params) const;
291
292 template<typename ValueType>
293 void writeBinaryGridImpl(const escript::Data& in,
294 const std::string& filename, int byteOrder) const;
295
296 dim_t findNode(const double *coords) const;
297
298
299 escript::Data randomFillWorker(const escript::DataTypes::ShapeType& shape,
300 long seed, const boost::python::tuple& filter) const;
301
302 #ifdef ESYS_MPI
303 /// corner neighbours' ranks
304 int neighbour_ranks[8];
305
306 /// corner neighbours' existence
307 bool neighbour_exists[8];
308 #endif
309
310 /// total number of elements in each dimension
311 dim_t m_gNE[3];
312
313 /// origin of domain
314 double m_origin[3];
315
316 /// side lengths of domain
317 double m_length[3];
318
319 /// grid spacings / cell sizes of domain
320 double m_dx[3];
321
322 /// number of spatial subdivisions
323 int m_NX[3];
324
325 /// number of elements for this rank in each dimension including shared
326 dim_t m_NE[3];
327
328 /// number of nodes for this rank in each dimension
329 dim_t m_NN[3];
330
331 /// first node on this rank is at (offset0,offset1,offset2) in global mesh
332 dim_t m_offset[3];
333
334 /// number of face elements per edge (left, right, bottom, top, front, back)
335 dim_t m_faceCount[6];
336
337
338 /// vector of sample reference identifiers
339 IndexVector m_dofId;
340 IndexVector m_nodeId;
341 IndexVector m_elementId;
342
343 // vector with first node id on each rank
344 IndexVector m_nodeDistribution;
345
346 #ifdef USE_RIPLEY
347 mutable RipleyCoupler *coupler;
348 #endif
349 };
350
351 ////////////////////////////// inline methods ////////////////////////////////
352 inline dim_t Brick::getDofOfNode(dim_t node) const
353 {
354 return m_nodeId[node];
355 }
356
357 inline dim_t Brick::getNumDataPointsGlobal() const
358 {
359 return (m_gNE[0]*m_order+1)*(m_gNE[1]*m_order+1)*(m_gNE[2]*m_order+1);
360 }
361
362 inline double Brick::getLocalCoordinate(index_t index, int dim) const
363 {
364 EsysAssert((dim>=0 && dim<m_numDim), "'dim' out of bounds");
365 EsysAssert((index>=0 && index<m_NN[dim]), "'index' out of bounds");
366 return m_origin[dim] //origin
367 + m_dx[dim]*(m_offset[dim] + index/m_order //elements
368 + point_locations[m_order-2][index%m_order]); //quads
369 }
370
371 inline boost::python::tuple Brick::getGridParameters() const
372 {
373 return boost::python::make_tuple(
374 boost::python::make_tuple(m_origin[0], m_origin[1], m_origin[2]),
375 boost::python::make_tuple(m_dx[0], m_dx[1], m_dx[2]),
376 boost::python::make_tuple(m_gNE[0], m_gNE[1], m_gNE[2]));
377 }
378
379 //protected
380 inline dim_t Brick::getNumDOF() const //global points
381 {
382 return getNumNodes();
383 }
384
385 //protected
386 inline dim_t Brick::getNumNodes() const //points per rank
387 {
388 return m_NN[0] * m_NN[1] * m_NN[2];
389 }
390
391 //protected
392 inline dim_t Brick::getNumElements() const
393 {
394 return m_NE[0]*m_NE[1]*m_NE[2];
395 }
396
397 } // end of namespace speckley
398
399 #endif // __Speckley_BRICK_H__
400

  ViewVC Help
Powered by ViewVC 1.1.26