libMesh
fe_xyz_map.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "libmesh/fe_xyz_map.h"
19 #include "libmesh/elem.h"
20 
21 namespace libMesh
22 {
23 
24 void FEXYZMap::compute_face_map(int dim, const std::vector<Real> & qw, const Elem * side)
25 {
26  libmesh_assert(side);
27 
28  LOG_SCOPE("compute_face_map()", "FEXYZMap");
29 
30  // The number of quadrature points.
31  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
32 
33  switch(dim)
34  {
35  case 2:
36  {
37 
38  // Resize the vectors to hold data at the quadrature points
39  {
40  this->xyz.resize(n_qp);
41  this->dxyzdxi_map.resize(n_qp);
42  this->d2xyzdxi2_map.resize(n_qp);
43  this->tangents.resize(n_qp);
44  this->normals.resize(n_qp);
45  this->curvatures.resize(n_qp);
46 
47  this->JxW.resize(n_qp);
48  }
49 
50  // Clear the entities that will be summed
51  // Compute the tangent & normal at the quadrature point
52  for (unsigned int p=0; p<n_qp; p++)
53  {
54  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
55  this->xyz[p].zero();
56  this->dxyzdxi_map[p].zero();
57  this->d2xyzdxi2_map[p].zero();
58  }
59 
60  // compute x, dxdxi at the quadrature points
61  for (std::size_t i=0; i<this->psi_map.size(); i++) // sum over the nodes
62  {
63  const Point & side_point = side->point(i);
64 
65  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
66  {
67  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
68  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
69  this->d2xyzdxi2_map[p].add_scaled(side_point, this->d2psidxi2_map[i][p]);
70  }
71  }
72 
73  // Compute the tangent & normal at the quadrature point
74  for (unsigned int p=0; p<n_qp; p++)
75  {
76  const Point n(this->dxyzdxi_map[p](1), -this->dxyzdxi_map[p](0), 0.);
77 
78  this->normals[p] = n.unit();
79  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
80 #if LIBMESH_DIM == 3 // Only good in 3D space
81  this->tangents[p][1] = this->dxyzdxi_map[p].cross(n).unit();
82 #endif
83  // The curvature is computed via the familiar Frenet formula:
84  // curvature = [d^2(x) / d (xi)^2] dot [normal]
85  // For a reference, see:
86  // F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill, p. 310
87  //
88  // Note: The sign convention here is different from the
89  // 3D case. Concave-upward curves (smiles) have a positive
90  // curvature. Concave-downward curves (frowns) have a
91  // negative curvature. Be sure to take that into account!
92  const Real numerator = this->d2xyzdxi2_map[p] * this->normals[p];
93  const Real denominator = this->dxyzdxi_map[p].norm_sq();
94  libmesh_assert_not_equal_to (denominator, 0);
95  this->curvatures[p] = numerator / denominator;
96  }
97 
98  // compute the jacobian at the quadrature points
99  for (unsigned int p=0; p<n_qp; p++)
100  {
101  const Real the_jac = std::sqrt(this->dxdxi_map(p)*this->dxdxi_map(p) +
102  this->dydxi_map(p)*this->dydxi_map(p));
103 
104  libmesh_assert_greater (the_jac, 0.);
105 
106  this->JxW[p] = the_jac*qw[p];
107  }
108 
109  break;
110  }
111 
112  case 3:
113  {
114  // Resize the vectors to hold data at the quadrature points
115  {
116  this->xyz.resize(n_qp);
117  this->dxyzdxi_map.resize(n_qp);
118  this->dxyzdeta_map.resize(n_qp);
119  this->d2xyzdxi2_map.resize(n_qp);
120  this->d2xyzdxideta_map.resize(n_qp);
121  this->d2xyzdeta2_map.resize(n_qp);
122  this->tangents.resize(n_qp);
123  this->normals.resize(n_qp);
124  this->curvatures.resize(n_qp);
125 
126  this->JxW.resize(n_qp);
127  }
128 
129  // Clear the entities that will be summed
130  for (unsigned int p=0; p<n_qp; p++)
131  {
132  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
133  this->xyz[p].zero();
134  this->dxyzdxi_map[p].zero();
135  this->dxyzdeta_map[p].zero();
136  this->d2xyzdxi2_map[p].zero();
137  this->d2xyzdxideta_map[p].zero();
138  this->d2xyzdeta2_map[p].zero();
139  }
140 
141  // compute x, dxdxi at the quadrature points
142  for (std::size_t i=0; i<this->psi_map.size(); i++) // sum over the nodes
143  {
144  const Point & side_point = side->point(i);
145 
146  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
147  {
148  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
149  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
150  this->dxyzdeta_map[p].add_scaled (side_point, this->dpsideta_map[i][p]);
151  this->d2xyzdxi2_map[p].add_scaled (side_point, this->d2psidxi2_map[i][p]);
152  this->d2xyzdxideta_map[p].add_scaled(side_point, this->d2psidxideta_map[i][p]);
153  this->d2xyzdeta2_map[p].add_scaled (side_point, this->d2psideta2_map[i][p]);
154  }
155  }
156 
157  // Compute the tangents, normal, and curvature at the quadrature point
158  for (unsigned int p=0; p<n_qp; p++)
159  {
160  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
161  this->normals[p] = n.unit();
162  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
163  this->tangents[p][1] = n.cross(this->dxyzdxi_map[p]).unit();
164 
165  // Compute curvature using the typical nomenclature
166  // of the first and second fundamental forms.
167  // For reference, see:
168  // 1) http://mathworld.wolfram.com/MeanCurvature.html
169  // (note -- they are using inward normal)
170  // 2) F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill
171  const Real L = -this->d2xyzdxi2_map[p] * this->normals[p];
172  const Real M = -this->d2xyzdxideta_map[p] * this->normals[p];
173  const Real N = -this->d2xyzdeta2_map[p] * this->normals[p];
174  const Real E = this->dxyzdxi_map[p].norm_sq();
175  const Real F = this->dxyzdxi_map[p] * this->dxyzdeta_map[p];
176  const Real G = this->dxyzdeta_map[p].norm_sq();
177 
178  const Real numerator = E*N -2.*F*M + G*L;
179  const Real denominator = E*G - F*F;
180  libmesh_assert_not_equal_to (denominator, 0.);
181  this->curvatures[p] = 0.5*numerator/denominator;
182  }
183 
184  // compute the jacobian at the quadrature points, see
185  // http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/thtoc.htm
186  for (unsigned int p=0; p<n_qp; p++)
187  {
188  const Real g11 = (this->dxdxi_map(p)*this->dxdxi_map(p) +
189  this->dydxi_map(p)*this->dydxi_map(p) +
190  this->dzdxi_map(p)*this->dzdxi_map(p));
191 
192  const Real g12 = (this->dxdxi_map(p)*this->dxdeta_map(p) +
193  this->dydxi_map(p)*this->dydeta_map(p) +
194  this->dzdxi_map(p)*this->dzdeta_map(p));
195 
196  const Real g21 = g12;
197 
198  const Real g22 = (this->dxdeta_map(p)*this->dxdeta_map(p) +
199  this->dydeta_map(p)*this->dydeta_map(p) +
200  this->dzdeta_map(p)*this->dzdeta_map(p));
201 
202 
203  const Real the_jac = std::sqrt(g11*g22 - g12*g21);
204 
205  libmesh_assert_greater (the_jac, 0.);
206 
207  this->JxW[p] = the_jac*qw[p];
208  }
209 
210  break;
211  }
212  default:
213  libmesh_error_msg("Invalid dim = " << dim);
214 
215  } // switch(dim)
216 }
217 
218 } // namespace libMesh
virtual void compute_face_map(int dim, const std::vector< Real > &qw, const Elem *side) libmesh_override
Special implementation for XYZ finite elements.
Definition: fe_xyz_map.C:24
std::vector< std::vector< Real > > dpsideta_map
Map for the derivative of the side function, d(psi)/d(eta).
Definition: fe_map.h:820
std::vector< std::vector< Real > > psi_map
Map for the side shape functions, psi.
Definition: fe_map.h:808
std::vector< std::vector< Real > > d2psideta2_map
Map for the second derivatives (in eta) of the side shape functions.
Definition: fe_map.h:841
unsigned int dim
Real dxdxi_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:546
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
The libMesh namespace provides an interface to certain functionality in the library.
Real dydxi_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:554
std::vector< std::vector< Real > > d2psidxideta_map
Map for the second (cross) derivatives in xi, eta of the side shape functions.
Definition: fe_map.h:834
std::vector< RealGradient > d2xyzdxideta_map
Vector of mixed second partial derivatives in xi-eta: d^2(x)/d(xi)d(eta) d^2(y)/d(xi)d(eta) d^2(z)/d(...
Definition: fe_map.h:645
libmesh_assert(j)
std::vector< RealGradient > dxyzdxi_map
Vector of partial derivatives: d(x)/d(xi), d(y)/d(xi), d(z)/d(xi)
Definition: fe_map.h:621
std::vector< RealGradient > d2xyzdeta2_map
Vector of second partial derivatives in eta: d^2(x)/d(eta)^2.
Definition: fe_map.h:651
std::vector< RealGradient > d2xyzdxi2_map
Vector of second partial derivatives in xi: d^2(x)/d(xi)^2, d^2(y)/d(xi)^2, d^2(z)/d(xi)^2.
Definition: fe_map.h:639
Real dzdxi_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:562
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Definition: type_vector.h:879
Real dxdeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:570
Real dydeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:578
std::vector< Real > curvatures
The mean curvature (= one half the sum of the principal curvatures) on the boundary at the quadrature...
Definition: fe_map.h:858
std::vector< std::vector< Real > > d2psidxi2_map
Map for the second derivatives (in xi) of the side shape functions.
Definition: fe_map.h:827
TypeVector< T > unit() const
Definition: type_vector.C:37
std::vector< Real > JxW
Jacobian*Weight values at quadrature points.
Definition: fe_map.h:868
std::vector< Point > normals
Normal vectors on boundary at quadrature points.
Definition: fe_map.h:851
std::vector< Point > xyz
The spatial locations of the quadrature points.
Definition: fe_map.h:615
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
std::vector< std::vector< Real > > dpsidxi_map
Map for the derivative of the side functions, d(psi)/d(xi).
Definition: fe_map.h:814
std::vector< std::vector< Point > > tangents
Tangent vectors on boundary at quadrature points.
Definition: fe_map.h:846
std::vector< RealGradient > dxyzdeta_map
Vector of partial derivatives: d(x)/d(eta), d(y)/d(eta), d(z)/d(eta)
Definition: fe_map.h:627
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
Real dzdeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:586