1 /**
2 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain rights in this software.
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 */
15
16 #ifndef SWEPT_ELEMENT_DATA_HPP
17 #define SWEPT_ELEMENT_DATA_HPP
18
19 //
20 // Class: SweptElementData
21 //
22 // Purpose: represent a rectangular element of mesh
23 //
24 // A SweptElement represents a rectangular element of mesh, including both vertices and
25 // elements, and the parametric space used to address that element. Vertex data,
26 // i.e. coordinates, may not be stored directly in the element, but the element returns
27 // information about the vertex handles of vertices in the element. Vertex and element
28 // handles associated with the element are each contiguous.
29
30 #include "SequenceData.hpp"
31 #include "moab/HomXform.hpp"
32 #include "moab/CN.hpp"
33 #include "SweptVertexData.hpp"
34 #include "Internals.hpp"
35 #include "moab/Range.hpp"
36
37 #include <vector>
38 #include <algorithm>
39
40 namespace moab {
41
42 class SweptElementData : public SequenceData
43 {
44
45 //! structure to hold references to bounding vertex blocks
46 class VertexDataRef
47 {
48 private:
49 HomCoord minmax[2];
50 HomXform xform, invXform;
51 SweptVertexData *srcSeq;
52 public:
53 friend class SweptElementData;
54
55 VertexDataRef(const HomCoord &min, const HomCoord &max,
56 const HomXform &tmp_xform, SweptVertexData *this_seq);
57
58 bool contains(const HomCoord &coords) const;
59 };
60
61 private:
62
63 //! parameter min/max/stride, in homogeneous coords ijkh
64 HomCoord elementParams[3];
65
66 //! difference between max and min params plus one (i.e. # VERTICES in
67 //! each parametric direction)
68 int dIJK[3];
69
70 //! difference between max and min params (i.e. # ELEMENTS in
71 //! each parametric direction)
72 int dIJKm1[3];
73
74 //! bare constructor, so compiler doesn't create one for me
75 SweptElementData();
76
77 //! list of bounding vertex blocks
78 std::vector<VertexDataRef> vertexSeqRefs;
79
80 public:
81
82 //! constructor
83 SweptElementData( EntityHandle start_handle,
84 const int imin, const int jmin, const int kmin,
85 const int imax, const int jmax, const int kmax,
86 const int* Cq );
87
88 virtual ~SweptElementData();
89
90 //! get handle of vertex at homogeneous coords
91 inline EntityHandle get_vertex(const HomCoord &coords) const;
get_vertex(int i,int j,int k) const92 inline EntityHandle get_vertex(int i, int j, int k) const
93 { return get_vertex(HomCoord(i,j,k)); }
94
95 //! get handle of element at i, j, k
96 EntityHandle get_element(const int i, const int j, const int k) const;
97
98 //! get min params for this element
99 const HomCoord &min_params() const;
100
101 //! get max params for this element
102 const HomCoord &max_params() const;
103
104 //! get the number of vertices in each direction, inclusive
105 void param_extents(int &di, int &dj, int &dk) const;
106
107 //! given a handle, get the corresponding parameters
108 ErrorCode get_params(const EntityHandle ehandle,
109 int &i, int &j, int &k) const;
110
111 //! convenience functions for parameter extents
i_min() const112 int i_min() const {return (elementParams[0].hom_coord())[0];}
j_min() const113 int j_min() const {return (elementParams[0].hom_coord())[1];}
k_min() const114 int k_min() const {return (elementParams[0].hom_coord())[2];}
i_max() const115 int i_max() const {return (elementParams[1].hom_coord())[0];}
j_max() const116 int j_max() const {return (elementParams[1].hom_coord())[1];}
k_max() const117 int k_max() const {return (elementParams[1].hom_coord())[2];}
118
119 //! test the bounding vertex sequences and determine whether they fully
120 //! define the vertices covering this element block's parameter space
121 bool boundary_complete() const;
122
123 //! test whether this sequence contains these parameters
124 inline bool contains(const HomCoord &coords) const;
125
126 //! get connectivity of an entity given entity's parameters
127 inline ErrorCode get_params_connectivity(const int i, const int j, const int k,
128 std::vector<EntityHandle>& connectivity) const;
129
130 //! add a vertex seq ref to this element sequence;
131 //! if bb_input is true, bounding box (in eseq-local coords) of vseq being added
132 //! is input in bb_min and bb_max (allows partial sharing of vseq rather than the whole
133 //! vseq); if it's false, the whole vseq is referenced and the eseq-local coordinates
134 //! is computed from the transformed bounding box of the vseq
135 ErrorCode add_vsequence(SweptVertexData *vseq,
136 const HomCoord &p1, const HomCoord &q1,
137 const HomCoord &p2, const HomCoord &q2,
138 const HomCoord &p3, const HomCoord &q3,
139 bool bb_input = false,
140 const HomCoord &bb_min = HomCoord::unitv[0],
141 const HomCoord &bb_max = HomCoord::unitv[0]);
142
143
144 SequenceData* subset( EntityHandle start,
145 EntityHandle end,
146 const int* sequence_data_sizes,
147 const int* tag_data_sizes ) const;
148
149 static EntityID calc_num_entities( EntityHandle start_handle,
150 int irange,
151 int jrange,
152 int krange );
153
154 unsigned long get_memory_use() const;
155 };
156
get_element(const int i,const int j,const int k) const157 inline EntityHandle SweptElementData::get_element(const int i, const int j, const int k) const
158 {
159 return start_handle() + (i-i_min()) + (j-j_min())*dIJKm1[0] + (k-k_min())*dIJKm1[0]*dIJKm1[1];
160 }
161
min_params() const162 inline const HomCoord &SweptElementData::min_params() const
163 {
164 return elementParams[0];
165 }
166
max_params() const167 inline const HomCoord &SweptElementData::max_params() const
168 {
169 return elementParams[1];
170 }
171
172 //! get the number of vertices in each direction, inclusive
param_extents(int & di,int & dj,int & dk) const173 inline void SweptElementData::param_extents(int &di, int &dj, int &dk) const
174 {
175 di = dIJK[0];
176 dj = dIJK[1];
177 dk = dIJK[2];
178 }
179
get_params(const EntityHandle ehandle,int & i,int & j,int & k) const180 inline ErrorCode SweptElementData::get_params(const EntityHandle ehandle,
181 int &i, int &j, int &k) const
182 {
183 if (TYPE_FROM_HANDLE(ehandle) != TYPE_FROM_HANDLE(start_handle())) return MB_FAILURE;
184
185 int hdiff = ehandle - start_handle();
186
187 // use double ?: test below because on some platforms, both sides of the : are
188 // evaluated, and if dIJKm1[1] is zero, that'll generate a divide-by-zero
189 k = (dIJKm1[1] > 0 ? hdiff / (dIJKm1[1] > 0 ? dIJKm1[0]*dIJKm1[1] : 1) : 0);
190 j = (hdiff - (k*dIJKm1[0]*dIJKm1[1])) / dIJKm1[0];
191 i = hdiff % dIJKm1[0];
192
193 k += elementParams[0].k();
194 j += elementParams[0].j();
195 i += elementParams[0].i();
196
197 return (ehandle >= start_handle() &&
198 ehandle < start_handle()+size() &&
199 i >= i_min() && i <= i_max() &&
200 j >= j_min() && j <= j_max() &&
201 k >= k_min() && k <= k_max()) ? MB_SUCCESS : MB_FAILURE;
202 }
203
contains(const HomCoord & temp) const204 inline bool SweptElementData::contains(const HomCoord &temp) const
205 {
206 // upper bound is < instead of <= because element params max is one less
207 // than vertex params max
208 return (temp >= elementParams[0] && temp < elementParams[1]);
209 }
210
contains(const HomCoord & coords) const211 inline bool SweptElementData::VertexDataRef::contains(const HomCoord &coords) const
212 {
213 return (minmax[0] <= coords && minmax[1] >= coords);
214 }
215
VertexDataRef(const HomCoord & this_min,const HomCoord & this_max,const HomXform & tmp_xform,SweptVertexData * this_seq)216 inline SweptElementData::VertexDataRef::VertexDataRef(const HomCoord &this_min, const HomCoord &this_max,
217 const HomXform &tmp_xform, SweptVertexData *this_seq)
218 : xform(tmp_xform), invXform(tmp_xform.inverse()), srcSeq(this_seq)
219 {
220 minmax[0] = HomCoord(this_min);
221 minmax[1] = HomCoord(this_max);
222 }
223
get_vertex(const HomCoord & coords) const224 inline EntityHandle SweptElementData::get_vertex(const HomCoord &coords) const
225 {
226 assert(boundary_complete());
227 for (std::vector<VertexDataRef>::const_iterator it = vertexSeqRefs.begin();
228 it != vertexSeqRefs.end(); ++it) {
229 if ((*it).minmax[0] <= coords && (*it).minmax[1] >= coords) {
230 // first get the vertex block-local parameters
231 HomCoord local_coords = coords / (*it).xform;
232
233 // now get the vertex handle for those coords
234 return (*it).srcSeq->get_vertex(local_coords);
235 }
236 }
237
238 // got here, it's an error
239 return 0;
240 }
241
add_vsequence(SweptVertexData * vseq,const HomCoord & p1,const HomCoord & q1,const HomCoord & p2,const HomCoord & q2,const HomCoord & p3,const HomCoord & q3,bool bb_input,const HomCoord & bb_min,const HomCoord & bb_max)242 inline ErrorCode SweptElementData::add_vsequence(SweptVertexData *vseq,
243 const HomCoord &p1, const HomCoord &q1,
244 const HomCoord &p2, const HomCoord &q2,
245 const HomCoord &p3, const HomCoord &q3,
246 bool bb_input,
247 const HomCoord &bb_min,
248 const HomCoord &bb_max)
249 {
250 // compute the transform given the vseq-local parameters and the mapping to
251 // this element sequence's parameters passed in minmax
252 HomXform M;
253 M.three_pt_xform(p1, q1, p2, q2, p3, q3);
254
255 // min and max in element seq's parameter system may not be same as those in
256 // vseq's system, so need to take min/max
257
258 HomCoord minmax[2];
259 if (bb_input) {
260 minmax[0] = bb_min;
261 minmax[1] = bb_max;
262 }
263 else {
264 minmax[0] = vseq->min_params() * M;
265 minmax[1] = vseq->max_params() * M;
266 }
267
268 // check against other vseq's to make sure they don't overlap
269 for (std::vector<VertexDataRef>::const_iterator vsit = vertexSeqRefs.begin();
270 vsit != vertexSeqRefs.end(); ++vsit)
271 if ((*vsit).contains(minmax[0]) || (*vsit).contains(minmax[1]))
272 return MB_FAILURE;
273
274 HomCoord tmp_min(std::min(minmax[0].i(), minmax[1].i()),
275 std::min(minmax[0].j(), minmax[1].j()),
276 std::min(minmax[0].k(), minmax[1].k()));
277 HomCoord tmp_max(std::max(minmax[0].i(), minmax[1].i()),
278 std::max(minmax[0].j(), minmax[1].j()),
279 std::max(minmax[0].k(), minmax[1].k()));
280
281
282 // set up a new vertex sequence reference
283 VertexDataRef tmp_seq_ref(tmp_min, tmp_max, M, vseq);
284
285 // add to the list
286 vertexSeqRefs.push_back(tmp_seq_ref);
287
288 return MB_SUCCESS;
289 }
290
get_params_connectivity(const int i,const int j,const int k,std::vector<EntityHandle> & connectivity) const291 inline ErrorCode SweptElementData::get_params_connectivity(const int i, const int j, const int k,
292 std::vector<EntityHandle>& connectivity) const
293 {
294 if (contains(HomCoord(i, j, k)) == false) return MB_FAILURE;
295
296 connectivity.push_back(get_vertex(i, j, k));
297 connectivity.push_back(get_vertex(i+1, j, k));
298 if (CN::Dimension(TYPE_FROM_HANDLE(start_handle())) < 2) return MB_SUCCESS;
299 connectivity.push_back(get_vertex(i+1, j+1, k));
300 connectivity.push_back(get_vertex(i, j+1, k));
301 if (CN::Dimension(TYPE_FROM_HANDLE(start_handle())) < 3) return MB_SUCCESS;
302 connectivity.push_back(get_vertex(i, j, k+1));
303 connectivity.push_back(get_vertex(i+1, j, k+1));
304 connectivity.push_back(get_vertex(i+1, j+1, k+1));
305 connectivity.push_back(get_vertex(i, j+1, k+1));
306 return MB_SUCCESS;
307 }
308
309 } // namespace moab
310
311 #endif
312