MFEM  v3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
eltrans.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 
13 #include <cmath>
14 #include "fem.hpp"
15 
16 namespace mfem
17 {
18 
20  JacobianIsEvaluated(0),
21  WeightIsEvaluated(0),
22  IntPoint(static_cast<IntegrationPoint *>(NULL)),
23  Attribute(-1),
24  ElementNo(-1)
25 {
26 
27 }
28 
30 {
31  switch (GeomType)
32  {
33  case Geometry::POINT : FElem = &PointFE; break;
34  case Geometry::SEGMENT : FElem = &SegmentFE; break;
35  case Geometry::TRIANGLE : FElem = &TriangleFE; break;
36  case Geometry::SQUARE : FElem = &QuadrilateralFE; break;
37  case Geometry::TETRAHEDRON : FElem = &TetrahedronFE; break;
38  case Geometry::CUBE : FElem = &HexahedronFE; break;
39  default:
40  MFEM_ABORT("unknown Geometry::Type!");
41  }
42  int dim = FElem->GetDim();
43  int dof = FElem->GetDof();
44  const IntegrationRule &nodes = FElem->GetNodes();
45  PointMat.SetSize(dim, dof);
46  for (int j = 0; j < dof; j++)
47  {
48  nodes.IntPoint(j).Get(&PointMat(0,j), dim);
49  }
50 }
51 
53 {
54  if (JacobianIsEvaluated) { return dFdx; }
55 
56  dshape.SetSize(FElem->GetDof(), FElem->GetDim());
57  dFdx.SetSize(PointMat.Height(), dshape.Width());
58 
59  FElem -> CalcDShape(*IntPoint, dshape);
60  Mult(PointMat, dshape, dFdx);
61 
63 
64  return dFdx;
65 }
66 
68 {
69  if (FElem->GetDim() == 0)
70  {
71  return 1.0;
72  }
74  {
75  return Wght;
76  }
77  Jacobian();
79  return (Wght = dFdx.Weight());
80 }
81 
83 {
84  switch (FElem->Space())
85  {
86  case FunctionSpace::Pk:
87  return (FElem->GetOrder()-1);
88  case FunctionSpace::Qk:
89  return (FElem->GetOrder());
90  default:
91  mfem_error("IsoparametricTransformation::OrderJ()");
92  }
93  return 0;
94 }
95 
97 {
98  switch (FElem->Space())
99  {
100  case FunctionSpace::Pk:
101  return (FElem->GetOrder() - 1) * FElem->GetDim();
102  case FunctionSpace::Qk:
103  return (FElem->GetOrder() * FElem->GetDim() - 1);
104  default:
105  mfem_error("IsoparametricTransformation::OrderW()");
106  }
107  return 0;
108 }
109 
111 {
112  if (FElem->Space() == fe->Space())
113  {
114  int k = FElem->GetOrder();
115  int d = FElem->GetDim();
116  int l = fe->GetOrder();
117  switch (fe->Space())
118  {
119  case FunctionSpace::Pk:
120  return ((k-1)*(d-1)+(l-1));
121  case FunctionSpace::Qk:
122  return (k*(d-1)+(l-1));
123  }
124  }
125  mfem_error("IsoparametricTransformation::OrderGrad(...)");
126  return 0;
127 }
128 
130  Vector &trans)
131 {
132  shape.SetSize(FElem->GetDof());
133  trans.SetSize(PointMat.Height());
134 
135  FElem -> CalcShape(ip, shape);
136  PointMat.Mult(shape, trans);
137 }
138 
140  DenseMatrix &tr)
141 {
142  int dof, n, dim, i, j, k;
143 
144  dim = PointMat.Height();
145  dof = FElem->GetDof();
146  n = ir.GetNPoints();
147 
148  shape.SetSize(dof);
149  tr.SetSize(dim, n);
150 
151  for (j = 0; j < n; j++)
152  {
153  FElem -> CalcShape (ir.IntPoint(j), shape);
154  for (i = 0; i < dim; i++)
155  {
156  tr(i, j) = 0.0;
157  for (k = 0; k < dof; k++)
158  {
159  tr(i, j) += PointMat(i, k) * shape(k);
160  }
161  }
162  }
163 }
164 
166  DenseMatrix &result)
167 {
168  result.SetSize(PointMat.Height(), matrix.Width());
169 
170  IntegrationPoint ip;
171  Vector col;
172 
173  for (int j = 0; j < matrix.Width(); j++)
174  {
175  ip.x = matrix(0, j);
176  if (matrix.Height() > 1)
177  {
178  ip.y = matrix(1, j);
179  if (matrix.Height() > 2)
180  {
181  ip.z = matrix(2, j);
182  }
183  }
184 
185  result.GetColumnReference(j, col);
186  Transform(ip, col);
187  }
188 }
189 
191  IntegrationPoint &ip)
192 {
193  const int max_iter = 16;
194  const double ref_tol = 1e-15;
195  const double phys_tol = 1e-15*pt.Normlinf();
196 
197  const int dim = FElem->GetDim();
198  const int sdim = PointMat.Height();
199  const int geom = FElem->GetGeomType();
200  IntegrationPoint xip, prev_xip;
201  double xd[3], yd[3], dxd[3], Jid[9];
202  Vector x(xd, dim), y(yd, sdim), dx(dxd, dim);
203  DenseMatrix Jinv(Jid, dim, sdim);
204  bool hit_bdr = false, prev_hit_bdr;
205 
206  // Use the center of the element as initial guess
207  xip = Geometries.GetCenter(geom);
208  xip.Get(xd, dim); // xip -> x
209 
210  for (int it = 0; it < max_iter; it++)
211  {
212  // Newton iteration: x := x + J(x)^{-1} [pt-F(x)]
213  // or when dim != sdim: x := x + [J^t.J]^{-1}.J^t [pt-F(x)]
214  Transform(xip, y);
215  subtract(pt, y, y); // y = pt-y
216  if (y.Normlinf() < phys_tol) { ip = xip; return 0; }
217  SetIntPoint(&xip);
218  CalcInverse(Jacobian(), Jinv);
219  Jinv.Mult(y, dx);
220  x += dx;
221  prev_xip = xip;
222  prev_hit_bdr = hit_bdr;
223  xip.Set(xd, dim); // x -> xip
224  // If xip is ouside project it on the boundary on the line segment
225  // between prev_xip and xip
226  hit_bdr = !Geometry::ProjectPoint(geom, prev_xip, xip);
227  if (dx.Normlinf() < ref_tol) { ip = xip; return 0; }
228  if (hit_bdr)
229  {
230  xip.Get(xd, dim); // xip -> x
231  if (prev_hit_bdr)
232  {
233  prev_xip.Get(dxd, dim); // prev_xip -> dx
234  subtract(x, dx, dx); // dx = xip - prev_xip
235  if (dx.Normlinf() < ref_tol) { return 1; }
236  }
237  }
238  }
239  ip = xip;
240  return 2;
241 }
242 
244  IntegrationPoint &ip2)
245 {
246  double vec[3];
247  Vector v (vec, Transf.GetPointMat().Height());
248 
249  Transf.Transform (ip1, v);
250  ip2.x = vec[0];
251  ip2.y = vec[1];
252  ip2.z = vec[2];
253 }
254 
256  IntegrationRule &ir2)
257 {
258  int i, n;
259 
260  n = ir1.GetNPoints();
261  for (i = 0; i < n; i++)
262  {
263  Transform (ir1.IntPoint(i), ir2.IntPoint(i));
264  }
265 }
266 
267 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:228
Abstract class for Finite Elements.
Definition: fe.hpp:44
void Set(const double *p, const int dim)
Definition: intrules.hpp:32
void Get(double *p, const int dim) const
Definition: intrules.hpp:45
int GetDim() const
Returns the space dimension for the finite element.
Definition: fe.hpp:82
Class for integration rule.
Definition: intrules.hpp:83
const IntegrationPoint * IntPoint
Definition: eltrans.hpp:28
void SetSize(int s)
Resize the vector if the new size is different.
Definition: vector.hpp:263
const Geometry::Type geom
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:451
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
Definition: operator.hpp:41
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:35
void SetIdentityTransformation(int GeomType)
Definition: eltrans.cpp:29
int GetOrder() const
Returns the order of the finite element.
Definition: fe.hpp:91
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
virtual int TransformBack(const Vector &, IntegrationPoint &)
Definition: eltrans.cpp:190
int Space() const
Returns the type of space on each element.
Definition: fe.hpp:94
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:649
virtual void Transform(const IntegrationPoint &, Vector &)
Definition: eltrans.cpp:129
const IntegrationPoint & GetCenter(int GeomType)
Definition: geom.hpp:57
double Weight() const
Definition: densemat.cpp:443
PointFiniteElement PointFE
Definition: point.cpp:30
Geometry Geometries
Definition: geom.cpp:544
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:231
int dim
Definition: ex3.cpp:47
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
Definition: operator.hpp:35
virtual int OrderGrad(const FiniteElement *fe)
order of adj(J)^t.grad(fi)
Definition: eltrans.cpp:110
TriLinear3DFiniteElement HexahedronFE
Definition: hexahedron.cpp:52
const IntegrationRule & GetNodes() const
Definition: fe.hpp:113
virtual const DenseMatrix & Jacobian()
Definition: eltrans.cpp:52
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:3012
int GetGeomType() const
Returns the geometry type:
Definition: fe.hpp:85
static bool ProjectPoint(int GeomType, const IntegrationPoint &beg, IntegrationPoint &end)
Definition: geom.cpp:317
int GetDof() const
Returns the degrees of freedom in the FE space.
Definition: fe.hpp:88
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:385
void mfem_error(const char *msg)
Definition: error.cpp:23
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:243
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:201
Linear3DFiniteElement TetrahedronFE
Linear2DFiniteElement TriangleFE
Definition: triangle.cpp:198
Class for integration point with weight.
Definition: intrules.hpp:25
IsoparametricTransformation Transf
Definition: eltrans.hpp:114
Vector data type.
Definition: vector.hpp:33
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:142
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:71
BiLinear2DFiniteElement QuadrilateralFE
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49