fe_xyz_map.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
20 namespace libMesh {
21 
22 void FEXYZMap::compute_face_map(int dim, const std::vector<Real>& qw, const Elem* side)
23 {
24  libmesh_assert(side);
25 
26  START_LOG("compute_face_map()", "FEXYZMap");
27 
28  // The number of quadrature points.
29  const unsigned int n_qp = libmesh_cast_int<unsigned int>(qw.size());
30 
31  switch(dim)
32  {
33  case 2:
34  {
35 
36  // Resize the vectors to hold data at the quadrature points
37  {
38  this->xyz.resize(n_qp);
39  this->dxyzdxi_map.resize(n_qp);
40  this->d2xyzdxi2_map.resize(n_qp);
41  this->tangents.resize(n_qp);
42  this->normals.resize(n_qp);
43  this->curvatures.resize(n_qp);
44 
45  this->JxW.resize(n_qp);
46  }
47 
48  // Clear the entities that will be summed
49  // Compute the tangent & normal at the quadrature point
50  for (unsigned int p=0; p<n_qp; p++)
51  {
52  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
53  this->xyz[p].zero();
54  this->dxyzdxi_map[p].zero();
55  this->d2xyzdxi2_map[p].zero();
56  }
57 
58  // compute x, dxdxi at the quadrature points
59  for (unsigned int i=0; i<this->psi_map.size(); i++) // sum over the nodes
60  {
61  const Point& side_point = side->point(i);
62 
63  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
64  {
65  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
66  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
67  this->d2xyzdxi2_map[p].add_scaled(side_point, this->d2psidxi2_map[i][p]);
68  }
69  }
70 
71  // Compute the tangent & normal at the quadrature point
72  for (unsigned int p=0; p<n_qp; p++)
73  {
74  const Point n(this->dxyzdxi_map[p](1), -this->dxyzdxi_map[p](0), 0.);
75 
76  this->normals[p] = n.unit();
77  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
78 #if LIBMESH_DIM == 3 // Only good in 3D space
79  this->tangents[p][1] = this->dxyzdxi_map[p].cross(n).unit();
80 #endif
81  // The curvature is computed via the familiar Frenet formula:
82  // curvature = [d^2(x) / d (xi)^2] dot [normal]
83  // For a reference, see:
84  // F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill, p. 310
85  //
86  // Note: The sign convention here is different from the
87  // 3D case. Concave-upward curves (smiles) have a positive
88  // curvature. Concave-downward curves (frowns) have a
89  // negative curvature. Be sure to take that into account!
90  const Real numerator = this->d2xyzdxi2_map[p] * this->normals[p];
91  const Real denominator = this->dxyzdxi_map[p].size_sq();
92  libmesh_assert_not_equal_to (denominator, 0);
93  this->curvatures[p] = numerator / denominator;
94  }
95 
96  // compute the jacobian at the quadrature points
97  for (unsigned int p=0; p<n_qp; p++)
98  {
99  const Real the_jac = std::sqrt(this->dxdxi_map(p)*this->dxdxi_map(p) +
100  this->dydxi_map(p)*this->dydxi_map(p));
101 
102  libmesh_assert_greater (the_jac, 0.);
103 
104  this->JxW[p] = the_jac*qw[p];
105  }
106 
107  break;
108  }
109 
110  case 3:
111  {
112  // Resize the vectors to hold data at the quadrature points
113  {
114  this->xyz.resize(n_qp);
115  this->dxyzdxi_map.resize(n_qp);
116  this->dxyzdeta_map.resize(n_qp);
117  this->d2xyzdxi2_map.resize(n_qp);
118  this->d2xyzdxideta_map.resize(n_qp);
119  this->d2xyzdeta2_map.resize(n_qp);
120  this->tangents.resize(n_qp);
121  this->normals.resize(n_qp);
122  this->curvatures.resize(n_qp);
123 
124  this->JxW.resize(n_qp);
125  }
126 
127  // Clear the entities that will be summed
128  for (unsigned int p=0; p<n_qp; p++)
129  {
130  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
131  this->xyz[p].zero();
132  this->dxyzdxi_map[p].zero();
133  this->dxyzdeta_map[p].zero();
134  this->d2xyzdxi2_map[p].zero();
135  this->d2xyzdxideta_map[p].zero();
136  this->d2xyzdeta2_map[p].zero();
137  }
138 
139  // compute x, dxdxi at the quadrature points
140  for (unsigned int i=0; i<this->psi_map.size(); i++) // sum over the nodes
141  {
142  const Point& side_point = side->point(i);
143 
144  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
145  {
146  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
147  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
148  this->dxyzdeta_map[p].add_scaled (side_point, this->dpsideta_map[i][p]);
149  this->d2xyzdxi2_map[p].add_scaled (side_point, this->d2psidxi2_map[i][p]);
150  this->d2xyzdxideta_map[p].add_scaled(side_point, this->d2psidxideta_map[i][p]);
151  this->d2xyzdeta2_map[p].add_scaled (side_point, this->d2psideta2_map[i][p]);
152  }
153  }
154 
155  // Compute the tangents, normal, and curvature at the quadrature point
156  for (unsigned int p=0; p<n_qp; p++)
157  {
158  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
159  this->normals[p] = n.unit();
160  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
161  this->tangents[p][1] = n.cross(this->dxyzdxi_map[p]).unit();
162 
163  // Compute curvature using the typical nomenclature
164  // of the first and second fundamental forms.
165  // For reference, see:
166  // 1) http://mathworld.wolfram.com/MeanCurvature.html
167  // (note -- they are using inward normal)
168  // 2) F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill
169  const Real L = -this->d2xyzdxi2_map[p] * this->normals[p];
170  const Real M = -this->d2xyzdxideta_map[p] * this->normals[p];
171  const Real N = -this->d2xyzdeta2_map[p] * this->normals[p];
172  const Real E = this->dxyzdxi_map[p].size_sq();
173  const Real F = this->dxyzdxi_map[p] * this->dxyzdeta_map[p];
174  const Real G = this->dxyzdeta_map[p].size_sq();
175 
176  const Real numerator = E*N -2.*F*M + G*L;
177  const Real denominator = E*G - F*F;
178  libmesh_assert_not_equal_to (denominator, 0.);
179  this->curvatures[p] = 0.5*numerator/denominator;
180  }
181 
182  // compute the jacobian at the quadrature points, see
183  // http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/thtoc.htm
184  for (unsigned int p=0; p<n_qp; p++)
185  {
186  const Real g11 = (this->dxdxi_map(p)*this->dxdxi_map(p) +
187  this->dydxi_map(p)*this->dydxi_map(p) +
188  this->dzdxi_map(p)*this->dzdxi_map(p));
189 
190  const Real g12 = (this->dxdxi_map(p)*this->dxdeta_map(p) +
191  this->dydxi_map(p)*this->dydeta_map(p) +
192  this->dzdxi_map(p)*this->dzdeta_map(p));
193 
194  const Real g21 = g12;
195 
196  const Real g22 = (this->dxdeta_map(p)*this->dxdeta_map(p) +
197  this->dydeta_map(p)*this->dydeta_map(p) +
198  this->dzdeta_map(p)*this->dzdeta_map(p));
199 
200 
201  const Real the_jac = std::sqrt(g11*g22 - g12*g21);
202 
203  libmesh_assert_greater (the_jac, 0.);
204 
205  this->JxW[p] = the_jac*qw[p];
206  }
207 
208  break;
209  }
210  default:
211  libmesh_error();
212 
213  } // switch(dim)
214 
215  STOP_LOG("compute_face_map()", "FEXYZMap");
216 
217  return;
218 }
219 
220 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo