1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
|
#ifndef FE_ELEMENT_H
#define FE_ELEMENT_H
#include <vector>
#include <string>
// ATC_Transfer headers
#include "MatrixLibrary.h"
#include "Array2D.h"
#include "ATC_TypeDefs.h"
namespace ATC {
enum ProjectionGuessType {
COORDINATE_ALIGNED=0,
CENTROID_LINEARIZED,
TWOD_ANALYTIC};
// Forward declarations
class FE_Interpolate;
/**
* @class FE_Element
* @brief Base class for a finite element holding info for canonical element
*/
class FE_Element {
public:
///////////////////////////////////////////////////////////////////////////
//
// CONSTRUCTOR AND DESTRUCTOR
FE_Element(const int nSD,
int numFaces,
int numNodes,
int numFaceNodes,
int numNodes1d);
virtual ~FE_Element();
///////////////////////////////////////////////////////////////////////////
//
// GETTERS
/** get number of spatial dimensions (almost always 3) */
int num_dims() { return nSD_; }
/** get number of element nodes */
int num_elt_nodes() { return numNodes_; }
/** get number of element nodes */
int num_elt_nodes_1d() { return numNodes1d_; }
/** get number of faces */
int num_faces() { return numFaces_; }
/** get number of face nodes */
int num_face_nodes() { return numFaceNodes_; }
// Getters for FE_Interpoate to have access to coordinates and connectivity
/** get canonical coordinates */
const DENS_MAT &local_coords() const { return localCoords_; }
/** get canonical coordinates in 1d */
DENS_VEC local_coords_1d() const;
/** get canonical connectivity of nodes and faces */
const Array2D<int> &local_face_conn() const { return localFaceConn_; }
/** return volume of the element */
double vol() const { return vol_; }
/** return area of a face */
double face_area() const { return faceArea_; }
// the following two are pass-throughs to the interpolate class, and
// can thus only be declared in the class body (or else the
// interpolate class is "incomplete" and cannot be referenced)
/** get number of integration points */
int num_ips() const;
/** get number of integration points */
int num_face_ips() const;
/** order of interpolation */
int order() const {return numNodes1d_;}
/** compute the quadrature for a given element type */
virtual void set_quadrature(FeIntQuadrature type) = 0;
/** return the set of 1d nodes that correspond to this node in 3d space */
void mapping(const int inode, std::vector<int> &mapping) const;
/** extract face coordinates from element coordinates */
void face_coordinates(const DENS_MAT &eltCoords,
const int faceID,
DENS_MAT &faceCoords) const;
/** set initial guess type for point in element search */
void set_projection_guess(ProjectionGuessType type)
{ projectionGuess_ = type;}
///////////////////////////////////////////////////////////////////////////
//
// GENERIC ELEMENT COMPUTATIONS
/** compute local coordinates from global */
virtual bool local_coordinates(const DENS_MAT &eltCoords,
const DENS_VEC &x,
DENS_VEC &xi) const;
/** location of local coordinates (0,0,0) */
virtual void centroid(const DENS_MAT &eltCoords,
DENS_VEC & centroid) const;
/** test if a specified element actually contains the given point */
virtual bool contains_point(const DENS_MAT &eltCoords,
const DENS_VEC &x) const;
/** check if element bounding box contains the given point */
bool range_check(const DENS_MAT &eltCoords, const DENS_VEC & x) const;
/** get the min and max coordinate of any point in an element in a
* dimension */
void bounds_in_dim(const DENS_MAT &eltCoords, const int dim,
double &min, double &max) const;
///////////////////////////////////////////////////////////////////////////
//
//PASS-THROUGHS TO INTERPOLATE CLASS
virtual void shape_function(const VECTOR & xi,
DENS_VEC &N) const;
/**
* compute shape functions at all ip's:
* indexed: N(ip,node)
* dN[nsd](ip,node)
* weights(ip)
*/
virtual void shape_function(const DENS_MAT eltCoords,
DENS_MAT &N,
std::vector<DENS_MAT> &dN,
DIAG_MAT &weights);
/**
* compute shape functions and derivatives at a single point,
* given the point and the element that contains it
* indexed: N(node)
*/
virtual void shape_function(const DENS_MAT eltCoords,
const VECTOR &x,
DENS_VEC &N);
/**
* compute shape functions and derivatives at a single point,
* given the point and the element that contains it
* indexed: N(node)
* dNdx(ip,nSD)
*/
virtual void shape_function(const DENS_MAT eltCoords,
const VECTOR &x,
DENS_VEC &N,
DENS_MAT &dNdx);
/**
* compute shape functions and derivatives at a single point,
* given the point and the element that contains it
* indexed:
* dNdx(ip,nSD)
*/
virtual void shape_function_derivatives(const DENS_MAT eltCoords,
const VECTOR &x,
DENS_MAT &dNdx);
/**
* compute shape functions at all face ip's:
* indexed: N(ip,node)
* n[nsd](ip,node)
* weights(ip)
*/
virtual void face_shape_function(const DENS_MAT &eltCoords,
const int faceID,
DENS_MAT &N,
DENS_MAT &n,
DIAG_MAT &weights);
/**
* compute shape functions at all face ip's:
* indexed: N(ip,node)
* dN[nsd](ip,node)
* Nn[nsd](ip,node)
* weights(ip)
*/
virtual void face_shape_function(const DENS_MAT &eltCoords,
const int faceID,
DENS_MAT &N,
std::vector<DENS_MAT> &dN,
std::vector<DENS_MAT> &Nn,
DIAG_MAT &weights);
/**
* compute normal vector from the specified face
* indexed: normal(nSD)
*/
virtual double face_normal(const DENS_MAT &eltCoords,
const int faceID,
int ip,
DENS_VEC &normal);
/**
* compute tangents to local coordinates
* indexed:
*/
virtual void tangents(const DENS_MAT &eltCoords,
const DENS_VEC &x,
std::vector<DENS_VEC> & tangents,
const bool normalize=false) const;
protected:
///////////////////////////////////////////////////////////////////////////
//
// HELPERS
/**
* generate the appropriate interpolation class
*/
FE_Interpolate *interpolate_factory(std::string interpolateType);
/** initial guess for local coordinates */
virtual void initial_local_coordinates(const DENS_MAT &eltCoords,
const DENS_VEC &x,
DENS_VEC &xiInitial) const;
///////////////////////////////////////////////////////////////////////////
//
// PROTECTED MEMBERS
// Currently used interpolation class
FE_Interpolate *feInterpolate_;
// Number of spatial dimensions
int nSD_;
// Number of faces, used for generic contains_point
int numFaces_;
// Number of element nodes
int numNodes_;
// Number of face nodes
int numFaceNodes_;
// Number of nodes in one dimension
int numNodes1d_;
// local coords of nodes: localCoords_(isd, ip)
DENS_MAT localCoords_;
// local face numbering
Array2D<int> localFaceConn_;
// volume of canonical element
double vol_;
// area of faces of canonical element
double faceArea_;
/** tolerance used in solving Newton's method for local coordinates */
double tolerance_;
ProjectionGuessType projectionGuess_;
};
/**
* @class FE_ElementHex
* @author Sean Laguna
* @brief 3D, linear 8-node hex element
*/
class FE_ElementHex : public FE_Element {
public:
FE_ElementHex(int numNodes,
int numFaceNodes,
int numNodes1d);
// Dump state info to disk for later restart (unimplemented)
void write_restart(FILE *);
~FE_ElementHex();
void set_quadrature(FeIntQuadrature type);
bool contains_point(const DENS_MAT &eltCoords,
const DENS_VEC &x) const;
};
/**
* @class FE_ElementRect
* @author Greg Wagner, amended by Sean Laguna
* @brief 3D, linear 8-node rectilinear hex element
*/
class FE_ElementRect : public FE_ElementHex {
public:
FE_ElementRect();
// Dump state info to disk for later restart (unimplemented)
void write_restart(FILE *);
~FE_ElementRect();
bool local_coordinates(const DENS_MAT &eltCoords,
const DENS_VEC &x,
DENS_VEC &xi) const;
protected:
virtual bool contains_point(const DENS_MAT &eltCoords,
const DENS_VEC &x) const;
};
/**
* @class FE_ElementTet
* @author Aaron Gable & Sean Laguna
* @brief 3D, linear 4-node tetrahedral element
*/
class FE_ElementTet : public FE_Element {
public:
FE_ElementTet(int numNodes,
int numFaceNodes,
int numNodes1d);
// Dump state info to disk for later restart (unimplemented)
void write_restart(FILE *);
~FE_ElementTet();
void set_quadrature(FeIntQuadrature type);
bool local_coordinates(const DENS_MAT &eltCoords,
const DENS_VEC &x,
DENS_VEC &xi) const;
bool contains_point(const DENS_MAT &eltCoords,
const DENS_VEC &x) const;
};
} // namespace ATC
#endif // FE_ELEMENT_H
|
