Branch data Line data Source code
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 SCD_ELEMENT_DATA_HPP
17 : : #define SCD_ELEMENT_DATA_HPP
18 : :
19 : : //
20 : : // Class: ScdElementData
21 : : //
22 : : // Purpose: represent a rectangular element of mesh
23 : : //
24 : : // A ScdElement 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 "ScdVertexData.hpp"
34 : : #include "Internals.hpp"
35 : : #include "moab/Range.hpp"
36 : :
37 : : #include <vector>
38 : : #include <algorithm>
39 : :
40 : : namespace moab
41 : : {
42 : :
43 : : class ScdElementData : public SequenceData
44 : : {
45 : :
46 : : //! structure to hold references to bounding vertex blocks
47 [ + + ]: 11325948 : class VertexDataRef
48 : : {
49 : : private:
50 : : HomCoord minmax[2];
51 : : HomXform xform, invXform;
52 : : ScdVertexData* srcSeq;
53 : :
54 : : public:
55 : : friend class ScdElementData;
56 : :
57 : : VertexDataRef( const HomCoord& min, const HomCoord& max, const HomXform& tmp_xform, ScdVertexData* this_seq );
58 : :
59 : : bool contains( const HomCoord& coords ) const;
60 : : };
61 : :
62 : : private:
63 : : //! parameter min/max/stride for vertices, in homogeneous coords ijkh; elements
64 : : //! are one less than max
65 : : HomCoord boxParams[3];
66 : :
67 : : //! difference between max and min params plus one (i.e. # VERTICES in
68 : : //! each parametric direction)
69 : : int dIJK[3];
70 : :
71 : : //! difference between max and min params (i.e. # ELEMENTS in
72 : : //! each parametric direction)
73 : : int dIJKm1[3];
74 : :
75 : : //! whether scd element rectangle is periodic in i and possibly j
76 : : int isPeriodic[2];
77 : :
78 : : //! bare constructor, so compiler doesn't create one for me
79 : : ScdElementData();
80 : :
81 : : //! list of bounding vertex blocks
82 : : std::vector< VertexDataRef > vertexSeqRefs;
83 : :
84 : : public:
85 : : //! constructor
86 : : ScdElementData( EntityHandle start_handle, const int imin, const int jmin, const int kmin, const int imax,
87 : : const int jmax, const int kmax, int* is_periodic );
88 : :
89 : : virtual ~ScdElementData();
90 : :
91 : : //! get handle of vertex at homogeneous coords
92 : : inline EntityHandle get_vertex( const HomCoord& coords ) const;
93 : 1409182 : inline EntityHandle get_vertex( int i, int j, int k ) const
94 : : {
95 [ + - ]: 1409182 : return get_vertex( HomCoord( i, j, k ) );
96 : : }
97 : :
98 : : //! get handle of element at i, j, k
99 : : EntityHandle get_element( const int i, const int j, const int k ) const;
100 : :
101 : : //! get min params for this element
102 : : const HomCoord& min_params() const;
103 : :
104 : : //! get max params for this element
105 : : const HomCoord& max_params() const;
106 : :
107 : : //! get the number of vertices in each direction, inclusive
108 : : void param_extents( int& di, int& dj, int& dk ) const;
109 : :
110 : : //! given a handle, get the corresponding parameters
111 : : ErrorCode get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const;
112 : :
113 : : //! return whether rectangle is periodic in i
114 : : inline int is_periodic_i() const
115 : : {
116 : : return isPeriodic[0];
117 : : };
118 : :
119 : : //! return whether rectangle is periodic in j
120 : : inline int is_periodic_j() const
121 : : {
122 : : return isPeriodic[1];
123 : : };
124 : :
125 : 13 : inline void is_periodic( int is_p[2] ) const
126 : : {
127 : 13 : is_p[0] = isPeriodic[0];
128 : 13 : is_p[1] = isPeriodic[1];
129 : 13 : }
130 : :
131 : : //! convenience functions for parameter extents
132 : 179690 : int i_min() const
133 : : {
134 : 179690 : return ( boxParams[0].hom_coord() )[0];
135 : : }
136 : 179690 : int j_min() const
137 : : {
138 : 179690 : return ( boxParams[0].hom_coord() )[1];
139 : : }
140 : 179690 : int k_min() const
141 : : {
142 : 179690 : return ( boxParams[0].hom_coord() )[2];
143 : : }
144 : 179690 : int i_max() const
145 : : {
146 : 179690 : return ( boxParams[1].hom_coord() )[0];
147 : : }
148 : 179690 : int j_max() const
149 : : {
150 : 179690 : return ( boxParams[1].hom_coord() )[1];
151 : : }
152 : 179690 : int k_max() const
153 : : {
154 : 179690 : return ( boxParams[1].hom_coord() )[2];
155 : : }
156 : :
157 : : //! test the bounding vertex sequences and determine whether they fully
158 : : //! define the vertices covering this element block's parameter space
159 : : bool boundary_complete() const;
160 : :
161 : : //! test whether this sequence contains these parameters
162 : : inline bool contains( const HomCoord& coords ) const;
163 : :
164 : : //! test whether *vertex parameterization* in this sequence contains these parameters
165 : : inline bool contains_vertex( const HomCoord& coords ) const;
166 : :
167 : : //! get connectivity of an entity given entity's parameters
168 : : inline ErrorCode get_params_connectivity( const int i, const int j, const int k,
169 : : std::vector< EntityHandle >& connectivity ) const;
170 : :
171 : : //! add a vertex seq ref to this element sequence;
172 : : //! if bb_input is true, bounding box (in eseq-local coords) of vseq being added
173 : : //! is input in bb_min and bb_max (allows partial sharing of vseq rather than the whole
174 : : //! vseq); if it's false, the whole vseq is referenced and the eseq-local coordinates
175 : : //! is computed from the transformed bounding box of the vseq
176 : : ErrorCode add_vsequence( ScdVertexData* vseq, const HomCoord& p1, const HomCoord& q1, const HomCoord& p2,
177 : : const HomCoord& q2, const HomCoord& p3, const HomCoord& q3, bool bb_input = false,
178 : : const HomCoord& bb_min = HomCoord::unitv[0], const HomCoord& bb_max = HomCoord::unitv[0] );
179 : :
180 : : SequenceData* subset( EntityHandle start, EntityHandle end, const int* sequence_data_sizes,
181 : : const int* tag_data_sizes ) const;
182 : :
183 : : static EntityID calc_num_entities( EntityHandle start_handle, int irange, int jrange, int krange,
184 : : int* is_periodic = NULL );
185 : :
186 : : unsigned long get_memory_use() const;
187 : : };
188 : :
189 : : inline EntityHandle ScdElementData::get_element( const int i, const int j, const int k ) const
190 : : {
191 : : return start_handle() + ( i - i_min() ) + ( j - j_min() ) * dIJKm1[0] + ( k - k_min() ) * dIJKm1[0] * dIJKm1[1];
192 : : }
193 : :
194 : 45 : inline const HomCoord& ScdElementData::min_params() const
195 : : {
196 : 45 : return boxParams[0];
197 : : }
198 : :
199 : 45 : inline const HomCoord& ScdElementData::max_params() const
200 : : {
201 : 45 : return boxParams[1];
202 : : }
203 : :
204 : : //! get the number of vertices in each direction, inclusive
205 : : inline void ScdElementData::param_extents( int& di, int& dj, int& dk ) const
206 : : {
207 : : di = dIJK[0];
208 : : dj = dIJK[1];
209 : : dk = dIJK[2];
210 : : }
211 : :
212 : 179690 : inline ErrorCode ScdElementData::get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const
213 : : {
214 [ - + ]: 179690 : if( TYPE_FROM_HANDLE( ehandle ) != TYPE_FROM_HANDLE( start_handle() ) ) return MB_FAILURE;
215 : :
216 : 179690 : int hdiff = ehandle - start_handle();
217 : :
218 : : // use double ?: test below because on some platforms, both sides of the : are
219 : : // evaluated, and if dIJKm1[1] is zero, that'll generate a divide-by-zero
220 [ + + ][ + - ]: 179690 : k = ( dIJKm1[1] > 0 ? hdiff / ( dIJKm1[1] > 0 ? dIJKm1[0] * dIJKm1[1] : 1 ) : 0 );
221 : 179690 : j = ( hdiff - ( k * dIJKm1[0] * dIJKm1[1] ) ) / dIJKm1[0];
222 : 179690 : i = hdiff % dIJKm1[0];
223 : :
224 : 179690 : k += boxParams[0].k();
225 : 179690 : j += boxParams[0].j();
226 : 179690 : i += boxParams[0].i();
227 : :
228 [ + - ][ + - ]: 539070 : return ( ehandle >= start_handle() && ehandle < start_handle() + size() && i >= i_min() && i <= i_max() &&
[ + - + - ]
229 [ + - ][ + - ]: 359380 : j >= j_min() && j <= j_max() && k >= k_min() && k <= k_max() )
[ + - ]
230 : : ? MB_SUCCESS
231 [ + - ]: 359380 : : MB_FAILURE;
232 : : }
233 : :
234 : 176184 : inline bool ScdElementData::contains( const HomCoord& temp ) const
235 : : {
236 : : // upper bound is < instead of <= because element params max is one less
237 : : // than vertex params max, except in case of 2d or 1d sequence
238 [ + - ][ + - ]: 176184 : return ( ( dIJKm1[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i() + dIJKm1[0] ) &&
[ + + ]
239 [ - + ][ + - ]: 176184 : ( ( !dIJKm1[1] && temp.j() == boxParams[1].j() ) ||
240 [ + - ][ + - ]: 528552 : ( dIJKm1[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j() + dIJKm1[1] ) ) &&
[ + - ][ + + ]
241 [ - + ][ + - ]: 176184 : ( ( !dIJKm1[2] && temp.k() == boxParams[1].k() ) ||
242 [ + - ][ + - ]: 352301 : ( dIJKm1[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k() + dIJKm1[2] ) ) );
243 : : }
244 : :
245 : : inline bool ScdElementData::contains_vertex( const HomCoord& temp ) const
246 : : {
247 : : // upper bound is < instead of <= because element params max is one less
248 : : // than vertex params max, except in case of 2d or 1d sequence
249 : : return ( ( dIJK[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i() + dIJK[0] ) &&
250 : : ( ( !dIJK[1] && temp.j() == boxParams[1].j() ) ||
251 : : ( dIJK[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j() + dIJK[1] ) ) &&
252 : : ( ( !dIJK[2] && temp.k() == boxParams[1].k() ) ||
253 : : ( dIJK[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k() + dIJK[2] ) ) );
254 : : }
255 : :
256 : 62435 : inline bool ScdElementData::VertexDataRef::contains( const HomCoord& coords ) const
257 : : {
258 [ + + ][ + + ]: 62435 : return ( minmax[0] <= coords && minmax[1] >= coords );
259 : : }
260 : :
261 : 39 : inline ScdElementData::VertexDataRef::VertexDataRef( const HomCoord& this_min, const HomCoord& this_max,
262 : : const HomXform& tmp_xform, ScdVertexData* this_seq )
263 [ + + ]: 117 : : xform( tmp_xform ), invXform( tmp_xform.inverse() ), srcSeq( this_seq )
264 : : {
265 [ + - ]: 39 : minmax[0] = HomCoord( this_min );
266 [ + - ]: 39 : minmax[1] = HomCoord( this_max );
267 : 39 : }
268 : :
269 : 1415702 : inline EntityHandle ScdElementData::get_vertex( const HomCoord& coords ) const
270 : : {
271 [ - + ]: 1415702 : assert( boundary_complete() );
272 [ + - ][ + - ]: 1421124 : for( std::vector< VertexDataRef >::const_iterator it = vertexSeqRefs.begin(); it != vertexSeqRefs.end(); ++it )
[ + - ]
273 : : {
274 [ + - ][ + - ]: 1421124 : if( ( *it ).minmax[0] <= coords && ( *it ).minmax[1] >= coords )
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ]
275 : : {
276 : : // first get the vertex block-local parameters
277 [ + - ][ + - ]: 1415702 : HomCoord local_coords = coords / ( *it ).xform;
278 : :
279 [ + - ][ + - ]: 1415702 : assert( ( *it ).srcSeq->contains( local_coords ) );
[ - + ]
280 : :
281 : : // now get the vertex handle for those coords
282 [ + - ][ + - ]: 1415702 : return ( *it ).srcSeq->get_vertex( local_coords );
283 : : }
284 : : }
285 : :
286 : : // got here, it's an error
287 : 1415702 : return 0;
288 : : }
289 : :
290 : 39 : inline ErrorCode ScdElementData::add_vsequence( ScdVertexData* vseq, const HomCoord& p1, const HomCoord& q1,
291 : : const HomCoord& p2, const HomCoord& q2, const HomCoord& p3,
292 : : const HomCoord& q3, bool bb_input, const HomCoord& bb_min,
293 : : const HomCoord& bb_max )
294 : : {
295 : : // compute the transform given the vseq-local parameters and the mapping to
296 : : // this element sequence's parameters passed in minmax
297 [ + - ]: 39 : HomXform M;
298 [ + - ]: 39 : M.three_pt_xform( p1, q1, p2, q2, p3, q3 );
299 : :
300 : : // min and max in element seq's parameter system may not be same as those in
301 : : // vseq's system, so need to take min/max
302 : :
303 [ + - ][ + + ]: 117 : HomCoord minmax[2];
304 [ + + ]: 39 : if( bb_input )
305 : : {
306 [ + - ]: 11 : minmax[0] = bb_min;
307 [ + - ]: 11 : minmax[1] = bb_max;
308 : : }
309 : : else
310 : : {
311 [ + - ][ + - ]: 28 : minmax[0] = vseq->min_params() * M;
[ + - ]
312 [ + - ][ + - ]: 28 : minmax[1] = vseq->max_params() * M;
[ + - ]
313 : : }
314 : :
315 : : // check against other vseq's to make sure they don't overlap
316 [ + - ][ + - ]: 53 : for( std::vector< VertexDataRef >::const_iterator vsit = vertexSeqRefs.begin(); vsit != vertexSeqRefs.end();
[ + - ][ + + ]
317 : : ++vsit )
318 [ + - ][ + - ]: 14 : if( ( *vsit ).contains( minmax[0] ) || ( *vsit ).contains( minmax[1] ) ) return MB_FAILURE;
[ + - ][ + - ]
[ + - ][ - + ]
[ - + ]
319 : :
320 [ + - ][ + - ]: 39 : HomCoord tmp_min( std::min( minmax[0].i(), minmax[1].i() ), std::min( minmax[0].j(), minmax[1].j() ),
[ + - ][ + - ]
[ + - ][ + - ]
321 [ + - ][ + - ]: 78 : std::min( minmax[0].k(), minmax[1].k() ) );
[ + - ][ + - ]
322 [ + - ][ + - ]: 39 : HomCoord tmp_max( std::max( minmax[0].i(), minmax[1].i() ), std::max( minmax[0].j(), minmax[1].j() ),
[ + - ][ + - ]
[ + - ][ + - ]
323 [ + - ][ + - ]: 78 : std::max( minmax[0].k(), minmax[1].k() ) );
[ + - ][ + - ]
324 : :
325 : : // set up a new vertex sequence reference
326 [ + - ]: 39 : VertexDataRef tmp_seq_ref( tmp_min, tmp_max, M, vseq );
327 : :
328 : : // add to the list
329 [ + - ]: 39 : vertexSeqRefs.push_back( tmp_seq_ref );
330 : :
331 : 39 : return MB_SUCCESS;
332 : : }
333 : :
334 : 176184 : inline ErrorCode ScdElementData::get_params_connectivity( const int i, const int j, const int k,
335 : : std::vector< EntityHandle >& connectivity ) const
336 : : {
337 [ + - ][ - + ]: 176184 : if( contains( HomCoord( i, j, k ) ) == false ) return MB_FAILURE;
338 : :
339 [ - + ]: 176184 : int ip1 = ( isPeriodic[0] ? ( i + 1 ) % dIJKm1[0] : i + 1 ),
340 [ - + ]: 176184 : jp1 = ( isPeriodic[1] ? ( j + 1 ) % dIJKm1[1] : j + 1 );
341 : :
342 [ + - ]: 176184 : connectivity.push_back( get_vertex( i, j, k ) );
343 [ + - ]: 176184 : connectivity.push_back( get_vertex( ip1, j, k ) );
344 [ + + ]: 176184 : if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 2 ) return MB_SUCCESS;
345 [ + - ]: 176173 : connectivity.push_back( get_vertex( ip1, jp1, k ) );
346 [ + - ]: 176173 : connectivity.push_back( get_vertex( i, jp1, k ) );
347 [ + + ]: 176173 : if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 3 ) return MB_SUCCESS;
348 [ + - ]: 176117 : connectivity.push_back( get_vertex( i, j, k + 1 ) );
349 [ + - ]: 176117 : connectivity.push_back( get_vertex( ip1, j, k + 1 ) );
350 [ + - ]: 176117 : connectivity.push_back( get_vertex( ip1, jp1, k + 1 ) );
351 [ + - ]: 176117 : connectivity.push_back( get_vertex( i, jp1, k + 1 ) );
352 : 176184 : return MB_SUCCESS;
353 : : }
354 : :
355 : : } // namespace moab
356 : :
357 : : #endif
|
