32 #ifndef CORE_BASIS_PRISMLINEARLGN_H
33 #define CORE_BASIS_PRISMLINEARLGN_H 1
104 const double dx = v1[0] - v0[0];
105 const double dy = v1[1] - v0[1];
106 const double dz = v1[2] - v0[2];
107 return sqrt(dx*dx+dy*dy+dz*dz);
123 template<
class VECTOR>
125 const unsigned div_per_unit,
126 VECTOR &coords)
const
128 typedef typename VECTOR::value_type VECTOR2;
130 coords.resize(div_per_unit + 1);
135 const double &p1x = v0[0];
136 const double &p1y = v0[1];
137 const double &p1z = v0[2];
138 const double dx = v1[0] - p1x;
139 const double dy = v1[1] - p1y;
140 const double dz = v1[2] - p1z;
142 for(
unsigned int i = 0; i <= div_per_unit; i++) {
143 const double d = (double)i / (
double)div_per_unit;
144 typename VECTOR::value_type &tmp = coords[i];
146 tmp[0] =
static_cast<typename VECTOR2::value_type
>(p1x + d * dx);
147 tmp[1] =
static_cast<typename VECTOR2::value_type
>(p1y + d * dy);
148 tmp[2] =
static_cast<typename VECTOR2::value_type
>(p1z + d * dz);
155 template<
class VECTOR>
157 const unsigned div_per_unit,
158 VECTOR &coords)
const
160 typedef typename VECTOR::value_type VECTOR2;
161 typedef typename VECTOR2::value_type VECTOR3;
164 coords.resize(fe * div_per_unit);
166 for(
unsigned int f = 0; f<fe; f++)
168 double *v0, *v1, *v2;
182 const double d = 1. / div_per_unit;
183 for(
unsigned int j = 0; j < div_per_unit; j++)
185 unsigned int k = j + f * div_per_unit;
186 const double dj = (double)j / (
double)div_per_unit;
188 coords[k].resize((div_per_unit - j) * 2 + 1);
189 typename VECTOR2::value_type &tmp = coords[k][e++];
191 tmp[0] =
static_cast<typename VECTOR3::value_type
>(v0[0] + dj * (v2[0] - v0[0]));
192 tmp[1] =
static_cast<typename VECTOR3::value_type
>(v0[1] + dj * (v2[1] - v0[1]));
193 tmp[2] =
static_cast<typename VECTOR3::value_type
>(v0[2] + dj * (v2[2] - v0[2]));
195 for(
unsigned int i = 0; i < div_per_unit - j; i++)
197 const double di = (double)i / (
double)div_per_unit;
198 typename VECTOR2::value_type &tmp1 = coords[j + f * div_per_unit][e++];
200 tmp1[0] =
static_cast<typename VECTOR3::value_type
>(v0[0] + (dj + d) * (v2[0] - v0[0]) + di * (v1[0] - v0[0]));
201 tmp1[1] =
static_cast<typename VECTOR3::value_type
>(v0[1] + (dj + d) * (v2[1] - v0[1]) + di * (v1[1] - v0[1]));
202 tmp1[2] =
static_cast<typename VECTOR3::value_type
>(v0[2] + (dj + d) * (v2[2] - v0[2]) + di * (v1[2] - v0[2]));
203 typename VECTOR2::value_type &tmp2 = coords[j + f * div_per_unit][e++];
205 tmp2[0] =
static_cast<typename VECTOR3::value_type
>(v0[0] + dj * (v2[0] - v0[0]) + (di + d) * (v1[0] - v0[0]));
206 tmp2[1] =
static_cast<typename VECTOR3::value_type
>(v0[1] + dj * (v2[1] - v0[1]) + (di + d) * (v1[1] - v0[1]));
207 tmp2[2] =
static_cast<typename VECTOR3::value_type
>(v0[2] + dj * (v2[2] - v0[2]) + (di + d) * (v1[2] - v0[2]));
216 template <
class ElemBasis>
219 typedef typename ElemBasis::value_type
T;
224 template <
class ElemData,
class VECTOR>
226 const T& value,
const ElemData &cd)
const
235 template <
class VECTOR>
249 template <
class ElemData,
class VECTOR>
251 VECTOR & guess)
const
253 double dist = DBL_MAX;
260 for (
int x = 1; x < end; x++)
262 coord[0] = x / (double) end;
263 for (
int y = 1; y < end; y++)
265 coord[1] = y / (double) end;
268 for (
int z = 1; z < end; z++)
270 coord[2] = z / (double) end;
275 pElem->derivate(coord, cd, derivs);
321 static T GaussianPoints[6][3];
322 static T GaussianWeights[6];
335 {1./6.,1./6., 0.211324865405}, {2./3.,1./6., 0.211324865405},
336 {1./6.,2./3., 0.211324865405}, {1./6.,1./6., 0.788675134595},
337 {2./3.,1./6., 0.788675134595}, {1./6.,2./3., 0.788675134595}};
340 T PrismGaussian2<T>::GaussianWeights[6] =
341 {1./6., 1./6., 1./6., 1./6., 1./6., 1./6.};
364 template<
class VECTOR>
366 { get_linear_weights(coords,w); }
368 template<
class VECTOR>
370 { get_linear_derivate_weights(coords,w); }
373 template <
class ElemData,
class VECTOR>
377 get_linear_weights(coords, w);
379 return (T)(w[0] * cd.node0() +
388 template <
class ElemData,
class VECTOR1,
class VECTOR2>
389 void derivate(
const VECTOR1 &coords,
const ElemData &cd,
390 VECTOR2 &derivs)
const
393 get_linear_derivate_weights(coords, w);
397 derivs[0] =
static_cast<typename VECTOR2::value_type
>(
405 derivs[1] =
static_cast<typename VECTOR2::value_type
>(
413 derivs[2] =
static_cast<typename VECTOR2::value_type
>(
419 w[17] * cd.node5()); }
422 template <
class ElemData,
class VECTOR>
424 const ElemData &cd)
const
427 return CL.
get_coords(
this, coords, value, cd);
431 template <
class ElemData>
434 return get_arc3d_length<CrvGaussian1<double> >(
this, edge, cd);
438 template <
class ElemData>
439 double get_area(
const unsigned face,
const ElemData &cd)
const
441 if (unit_faces[face][3]==-1)
442 return get_area3<TriGaussian2<double> >(
this, face, cd);
444 return get_area3<QuadGaussian2<double> >(
this, face, cd);
448 template <
class ElemData>
455 static const std::string type_name(
int n = -1);
469 ASSERT((n >= -1) && n <= 1);
477 static const std::string nm(
"PrismLinearLgn");
497 std::string(__FILE__),
void initial_guess(const ElemBasis *pElem, const T &val, const ElemData &cd, VECTOR &guess) const
find a reasonable initial guess
Definition: PrismLinearLgn.h:250
void approx_face(const unsigned face, const unsigned div_per_unit, VECTOR &coords) const
Definition: PrismLinearLgn.h:156
T value_type
Definition: PrismLinearLgn.h:357
Class for creating geometrical approximations of Prism meshes.
Definition: PrismLinearLgn.h:114
static int GaussianNum
Definition: PrismLinearLgn.h:320
double check_zero(const VECTOR &derivs, double epsilon=1e-7)
Definition: Locate.h:709
PrismLinearLgnUnitElement()
Definition: PrismLinearLgn.h:65
Class for describing unit geometry of PrismLinearLgn.
Definition: PrismLinearLgn.h:50
static SCISHARE int unit_edges[9][2]
References to vertices of unit edge.
Definition: PrismLinearLgn.h:55
static int vertices_of_face()
return number of vertices per face
Definition: PrismLinearLgn.h:93
static SCISHARE double unit_center[3]
Center of the unit element.
Definition: PrismLinearLgn.h:63
static SCISHARE double unit_vertices[6][3]
Parametric coordinates of vertices of unit edge.
Definition: PrismLinearLgn.h:53
Definition: Persistent.h:89
Definition: PrismLinearLgn.h:217
bool get_coords(const ElemBasis *pEB, VECTOR &coords, const T &value, const ElemData &cd) const
Definition: PrismLinearLgn.h:225
Class with weights and coordinates for 2nd order Gaussian integration.
Definition: PrismLinearLgn.h:317
virtual ~PrismLinearLgnUnitElement()
Definition: PrismLinearLgn.h:66
void get_derivate_weights(const VECTOR &coords, double *w) const
Definition: PrismLinearLgn.h:369
Definition: TypeDescription.h:45
#define SCISHARE
Definition: share.h:39
std::vector< const TypeDescription * > td_vec
Definition: TypeDescription.h:56
static const std::string type_name(int n=-1)
Definition: PrismLinearLgn.h:467
static int faces_of_cell()
return number of faces per cell
Definition: PrismLinearLgn.h:97
static T GaussianPoints[1][3]
Definition: PrismLinearLgn.h:295
ElemBasis::value_type T
Definition: Locate.h:549
virtual void io(Piostream &str)
Definition: PrismLinearLgn.h:507
#define ASSERT(condition)
Definition: Assert.h:110
Class for describing interfaces to basis elements.
Definition: Basis.h:48
virtual int begin_class(const std::string &name, int current_version)
Definition: Persistent.cc:143
const string find_type_name(float *)
Definition: TypeName.cc:63
static SCISHARE int unit_faces[5][4]
References to vertices of unit face.
Definition: PrismLinearLgn.h:57
static T GaussianWeights[1]
Definition: PrismLinearLgn.h:296
Class with weights and coordinates for 2nd order Gaussian integration.
Definition: PrismLinearLgn.h:291
static int number_of_vertices()
return number of vertices
Definition: PrismLinearLgn.h:77
const char * name[]
Definition: BoostGraphExampleTests.cc:87
void approx_edge(const unsigned edge, const unsigned div_per_unit, VECTOR &coords) const
Definition: PrismLinearLgn.h:124
Definition: StackVector.h:50
double get_arc_length(const unsigned edge, const ElemData &cd) const
get arc length for edge
Definition: PrismLinearLgn.h:432
static int GaussianNum
Definition: PrismLinearLgn.h:294
static double volume()
return volume
Definition: PrismLinearLgn.h:110
bool check_coords(const VECTOR &x) const
Definition: PrismLinearLgn.h:236
T interpolate(const VECTOR &coords, const ElemData &cd) const
get value at parametric coordinate
Definition: PrismLinearLgn.h:374
double get_area(const unsigned face, const ElemData &cd) const
get area
Definition: PrismLinearLgn.h:439
static double length(int edge)
Definition: PrismLinearLgn.h:100
static int number_of_edges()
return number of edges
Definition: PrismLinearLgn.h:85
void derivate(const VECTOR1 &coords, const ElemData &cd, VECTOR2 &derivs) const
get first derivative at parametric coordinate
Definition: PrismLinearLgn.h:389
Definition: PrismElementWeights.h:36
PrismApprox()
Definition: PrismLinearLgn.h:117
virtual ~PrismApprox()
Definition: PrismLinearLgn.h:118
virtual ~PrismLinearLgn()
Definition: PrismLinearLgn.h:360
Definition: PrismSamplingSchemes.h:42
static const std::string make_template_id(const std::string &templateName, const std::string &templateParam)
Definition: TypeName.h:62
static SCISHARE double unit_face_normals[5][3]
Normals of unit face.
Definition: PrismLinearLgn.h:59
static SCISHARE double unit_face_areas[5]
Precalculated area of faces.
Definition: PrismLinearLgn.h:61
static int polynomial_order()
Definition: PrismLinearLgn.h:362
static int dofs()
return degrees of freedom
Definition: PrismLinearLgn.h:89
PrismLinearLgn()
Definition: PrismLinearLgn.h:359
bool compare_distance(const T &interp, const T &val, double &cur_d, double dist)
Definition: Locate.h:667
virtual void end_class()
Definition: Persistent.cc:178
virtual ~PrismLocate()
Definition: PrismLinearLgn.h:222
static double domain_size()
return size of the domain
Definition: PrismLinearLgn.h:73
void resize(size_t size, const value_type &val=value_type())
Definition: StackVector.h:61
Class for handling of element of type prism with linear lagrangian interpolation. ...
Definition: PrismLinearLgn.h:348
PrismLocate()
Definition: PrismLinearLgn.h:221
ElemBasis::value_type T
Definition: PrismLinearLgn.h:219
bool get_coords(VECTOR &coords, const T &value, const ElemData &cd) const
get parametric coordinate for value within the element
Definition: PrismLinearLgn.h:423
const int PRISMLINEARLGN_VERSION
Definition: PrismLinearLgn.h:504
int n
Definition: eab.py:9
bool get_iterative(const ElemBasis *pEB, VECTOR &x, const T &value, const ElemData &cd) const
find value in interpolation for given value
Definition: Locate.h:559
double get_volume3(const ElemBasis *pEB, const ElemData &cd)
Definition: Locate.h:179
static int number_of_mesh_vertices()
return number of vertices in mesh
Definition: PrismLinearLgn.h:81
double get_volume(const ElemData &cd) const
get volume
Definition: PrismLinearLgn.h:449
static double area(int face)
return area
Definition: PrismLinearLgn.h:109
static int domain_dimension()
return dimension of domain
Definition: PrismLinearLgn.h:69
Definition: TypeDescription.h:49
void get_weights(const VECTOR &coords, double *w) const
Definition: PrismLinearLgn.h:365
const TypeDescription * get_type_description(Core::Basis::ConstantBasis< T > *)
Definition: Constant.h:209