MFEM  v4.6.0
Finite element discretization library
fe_ser.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 // Serendipity Finite Element classes
13 
14 #include "fe_ser.hpp"
15 #include "fe_fixed_order.hpp"
16 
17 namespace mfem
18 {
19 
20 using namespace std;
21 
23  : ScalarFiniteElement(2, Geometry::SQUARE, (p*p + 3*p +6) / 2, p,
24  FunctionSpace::Qk)
25 {
26  // Store the dof_map of the associated TensorBasisElement, which will be used
27  // to create the serendipity dof map. Its size is larger than the size of
28  // the serendipity element.
29  TensorBasisElement tbeTemp =
31  TensorBasisElement::DofMapType::Sr_DOF_MAP);
32  const Array<int> tp_dof_map = tbeTemp.GetDofMap();
33 
34  const double *cp = poly1d.ClosedPoints(p, BasisType::GaussLobatto);
35 
36  // Fixing the Nodes is exactly the same as the H1_QuadrilateralElement
37  // constructor except we only use those values of the associated tensor
38  // product dof_map that are <= the number of serendipity Dofs e.g. only DoFs
39  // 0-7 out of the 9 tensor product dofs (at quadratic order)
40  int o = 0;
41 
42  for (int j = 0; j <= p; j++)
43  {
44  for (int i = 0; i <= p; i++)
45  {
46  if (tp_dof_map[o] < Nodes.Size())
47  {
48  Nodes.IntPoint(tp_dof_map[o]).x = cp[i];
49  Nodes.IntPoint(tp_dof_map[o]).y = cp[j];
50  }
51  o++;
52  }
53  }
54 }
55 
57  Vector &shape) const
58 {
59  int p = (this)->GetOrder();
60  double x = ip.x, y = ip.y;
61 
63  Vector nodalX(p+1);
64  Vector nodalY(p+1);
65 
66  edgeNodalBasis.Eval(x, nodalX);
67  edgeNodalBasis.Eval(y, nodalY);
68 
69  // First, fix edge-based shape functions. Use a nodal interpolant for edge
70  // points, weighted by the linear function that vanishes on opposite edge.
71  for (int i = 0; i < p-1; i++)
72  {
73  shape(4 + 0*(p-1) + i) = (nodalX(i+1))*(1.-y); // south edge 0->1
74  shape(4 + 1*(p-1) + i) = (nodalY(i+1))*x; // east edge 1->2
75  shape(4 + 3*(p-1) - i - 1) = (nodalX(i+1)) * y; // north edge 3->2
76  shape(4 + 4*(p-1) - i - 1) = (nodalY(i+1)) * (1. - x); // west edge 0->3
77  }
78 
80  Vector bilinearsAtIP(4);
81  bilinear.CalcShape(ip, bilinearsAtIP);
82 
83  const double *edgePts(poly1d.ClosedPoints(p, BasisType::GaussLobatto));
84 
85  // Next, set the shape function associated with vertex V, evaluated at (x,y)
86  // to be: bilinear function associated to V, evaluated at (x,y) - sum (shape
87  // function at edge point P, weighted by bilinear function for V evaluated at
88  // P) where the sum is taken only for points P on edges incident to V.
89 
90  double vtx0fix =0;
91  double vtx1fix =0;
92  double vtx2fix =0;
93  double vtx3fix =0;
94  for (int i = 0; i<p-1; i++)
95  {
96  vtx0fix += (1-edgePts[i+1])*(shape(4 + i) +
97  shape(4 + 4*(p-1) - i - 1)); // bot+left edge
98  vtx1fix += (1-edgePts[i+1])*(shape(4 + 1*(p-1) + i) +
99  shape(4 + (p-2)-i)); // right+bot edge
100  vtx2fix += (1-edgePts[i+1])*(shape(4 + 2*(p-1) + i) +
101  shape(1 + 2*p-i)); // top+right edge
102  vtx3fix += (1-edgePts[i+1])*(shape(4 + 3*(p-1) + i) +
103  shape(3*p - i)); // left+top edge
104  }
105  shape(0) = bilinearsAtIP(0) - vtx0fix;
106  shape(1) = bilinearsAtIP(1) - vtx1fix;
107  shape(2) = bilinearsAtIP(2) - vtx2fix;
108  shape(3) = bilinearsAtIP(3) - vtx3fix;
109 
110  // Interior basis functions appear starting at order p=4. These are non-nodal
111  // bubble functions.
112  if (p > 3)
113  {
114  double *legX = new double[p-1];
115  double *legY = new double[p-1];
116  Poly_1D *storeLegendre = new Poly_1D();
117 
118  storeLegendre->CalcLegendre(p-2, x, legX);
119  storeLegendre->CalcLegendre(p-2, y, legY);
120 
121  int interior_total = 0;
122  for (int j = 4; j < p + 1; j++)
123  {
124  for (int k = 0; k < j-3; k++)
125  {
126  shape(4 + 4*(p-1) + interior_total)
127  = legX[k] * legY[j-4-k] * x * (1. - x) * y * (1. - y);
128  interior_total++;
129  }
130  }
131 
132  delete[] legX;
133  delete[] legY;
134  delete storeLegendre;
135  }
136 }
137 
139  DenseMatrix &dshape) const
140 {
141  int p = (this)->GetOrder();
142  double x = ip.x, y = ip.y;
143 
145  Vector nodalX(p+1);
146  Vector DnodalX(p+1);
147  Vector nodalY(p+1);
148  Vector DnodalY(p+1);
149 
150  edgeNodalBasis.Eval(x, nodalX, DnodalX);
151  edgeNodalBasis.Eval(y, nodalY, DnodalY);
152 
153  for (int i = 0; i < p-1; i++)
154  {
155  dshape(4 + 0*(p-1) + i,0) = DnodalX(i+1) * (1.-y);
156  dshape(4 + 0*(p-1) + i,1) = -nodalX(i+1);
157  dshape(4 + 1*(p-1) + i,0) = nodalY(i+1);
158  dshape(4 + 1*(p-1) + i,1) = DnodalY(i+1)*x;
159  dshape(4 + 3*(p-1) - i - 1,0) = DnodalX(i+1)*y;
160  dshape(4 + 3*(p-1) - i - 1,1) = nodalX(i+1);
161  dshape(4 + 4*(p-1) - i - 1,0) = -nodalY(i+1);
162  dshape(4 + 4*(p-1) - i - 1,1) = DnodalY(i+1) * (1.-x);
163  }
164 
166  DenseMatrix DbilinearsAtIP(4);
167  bilinear.CalcDShape(ip, DbilinearsAtIP);
168 
169  const double *edgePts(poly1d.ClosedPoints(p, BasisType::GaussLobatto));
170 
171  dshape(0,0) = DbilinearsAtIP(0,0);
172  dshape(0,1) = DbilinearsAtIP(0,1);
173  dshape(1,0) = DbilinearsAtIP(1,0);
174  dshape(1,1) = DbilinearsAtIP(1,1);
175  dshape(2,0) = DbilinearsAtIP(2,0);
176  dshape(2,1) = DbilinearsAtIP(2,1);
177  dshape(3,0) = DbilinearsAtIP(3,0);
178  dshape(3,1) = DbilinearsAtIP(3,1);
179 
180  for (int i = 0; i<p-1; i++)
181  {
182  dshape(0,0) -= (1-edgePts[i+1])*(dshape(4 + 0*(p-1) + i, 0) +
183  dshape(4 + 4*(p-1) - i - 1,0));
184  dshape(0,1) -= (1-edgePts[i+1])*(dshape(4 + 0*(p-1) + i, 1) +
185  dshape(4 + 4*(p-1) - i - 1,1));
186  dshape(1,0) -= (1-edgePts[i+1])*(dshape(4 + 1*(p-1) + i, 0) +
187  dshape(4 + (p-2)-i, 0));
188  dshape(1,1) -= (1-edgePts[i+1])*(dshape(4 + 1*(p-1) + i, 1) +
189  dshape(4 + (p-2)-i, 1));
190  dshape(2,0) -= (1-edgePts[i+1])*(dshape(4 + 2*(p-1) + i, 0) +
191  dshape(1 + 2*p-i, 0));
192  dshape(2,1) -= (1-edgePts[i+1])*(dshape(4 + 2*(p-1) + i, 1) +
193  dshape(1 + 2*p-i, 1));
194  dshape(3,0) -= (1-edgePts[i+1])*(dshape(4 + 3*(p-1) + i, 0) +
195  dshape(3*p - i, 0));
196  dshape(3,1) -= (1-edgePts[i+1])*(dshape(4 + 3*(p-1) + i, 1) +
197  dshape(3*p - i, 1));
198  }
199 
200  if (p > 3)
201  {
202  double *legX = new double[p-1];
203  double *legY = new double[p-1];
204  double *DlegX = new double[p-1];
205  double *DlegY = new double[p-1];
206  Poly_1D *storeLegendre = new Poly_1D();
207 
208  storeLegendre->CalcLegendre(p-2, x, legX, DlegX);
209  storeLegendre->CalcLegendre(p-2, y, legY, DlegY);
210 
211  int interior_total = 0;
212  for (int j = 4; j < p + 1; j++)
213  {
214  for (int k = 0; k < j-3; k++)
215  {
216  dshape(4 + 4*(p-1) + interior_total, 0) =
217  legY[j-4-k]*y*(1-y) * (DlegX[k]*x*(1-x) + legX[k]*(1-2*x));
218  dshape(4 + 4*(p-1) + interior_total, 1) =
219  legX[k]*x*(1-x) * (DlegY[j-4-k]*y*(1-y) + legY[j-4-k]*(1-2*y));
220  interior_total++;
221  }
222  }
223  delete[] legX;
224  delete[] legY;
225  delete[] DlegX;
226  delete[] DlegY;
227  delete storeLegendre;
228  }
229 }
230 
232  &Trans,
233  DenseMatrix &I) const
234 {
235  // For p<=4, the basis is nodal; for p>4, the quad-interior functions are
236  // non-nodal.
237  if (order <= 4)
238  {
239  NodalLocalInterpolation(Trans, I, *this);
240  }
241  else
242  {
243  ScalarLocalInterpolation(Trans, I, *this);
244  }
245 }
246 
247 }
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
STL namespace.
static void CalcLegendre(const int p, const double x, double *u)
Definition: fe_base.cpp:2085
const double * ClosedPoints(const int p, const int btype=BasisType::GaussLobatto)
Get coordinates of a closed (GaussLegendre) set of points if degree p.
Definition: fe_base.hpp:1071
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
const Array< int > & GetDofMap() const
Get an Array<int> that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe_base.hpp:1210
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Class for finite elements with basis functions that return scalar values.
Definition: fe_base.hpp:649
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
Class for computing 1D special polynomials and their associated basis functions.
Definition: fe_base.hpp:963
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition: fe_ser.cpp:56
Basis & GetBasis(const int p, const int btype)
Get a Poly_1D::Basis object of the given degree and BasisType, btype.
Definition: fe_base.cpp:2208
IntegrationRule Nodes
Definition: fe_base.hpp:246
void ScalarLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get matrix I "Interpolation" defined through local L2-projection in the space defined by the fine_fe...
Definition: fe_base.cpp:548
Poly_1D poly1d
Definition: fe.cpp:28
A 2D bi-linear element on a square with nodes at the vertices of the square.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_ser.cpp:138
Class for integration point with weight.
Definition: intrules.hpp:31
Class for evaluating 1D nodal, positive (Bernstein), or integrated (Gerritsma) bases.
Definition: fe_base.hpp:978
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
H1Ser_QuadrilateralElement(const int p)
Construct the H1Ser_QuadrilateralElement of order p.
Definition: fe_ser.cpp:22
Vector data type.
Definition: vector.hpp:58
Describes the function space on each element.
Definition: fe_base.hpp:216
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Definition: fe_ser.cpp:231
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
int order
Order/degree of the shape functions.
Definition: fe_base.hpp:243
void NodalLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get the matrix I that defines nodal interpolation between this element and the refined element fine_f...
Definition: fe_base.cpp:509