1 // Copyright (c) 2010-2021, 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 #include "mfem.hpp"
13 #include "unit_tests.hpp"
14
15 using namespace mfem;
16
testQuadratureInterpolator(const int dim,const int p,const int qpts,const QVectorLayout q_layout,const int nx,const int ny,const int nz)17 static bool testQuadratureInterpolator(const int dim,
18 const int p,
19 const int qpts,
20 const QVectorLayout q_layout,
21 const int nx, const int ny, const int nz)
22 {
23 // Keep for debugging purposes:
24 const bool verbose = false;
25 if (verbose)
26 {
27 std::cout << "testQuadratureInterpolator(dim=" << dim
28 << ",p=" << p
29 << ",q=" << qpts
30 << ",l=" << (q_layout == QVectorLayout::byNODES ?
31 "by_nodes" : "by_vdim")
32 << ",nx=" << nx
33 << ",ny=" << ny
34 << ",nz=" << nz
35 << ")" << std::endl;
36 }
37
38 const int vdim = dim;
39 const int seed = 0x100001b3;
40 const int ordering = Ordering::byNODES;
41
42 REQUIRE((dim == 2 || dim == 3));
43 Mesh mesh = dim == 1 ? Mesh::MakeCartesian1D(nx, Element::SEGMENT) :
44 dim == 2 ? Mesh::MakeCartesian2D(nx,ny, Element::QUADRILATERAL):
45 Mesh::MakeCartesian3D(nx,nx,nz, Element::HEXAHEDRON);
46
47 const H1_FECollection fec(p, dim);
48 FiniteElementSpace sfes(&mesh, &fec, 1, ordering);
49 FiniteElementSpace vfes(&mesh, &fec, vdim, ordering);
50
51 GridFunction x(&sfes);
52 x.Randomize(seed);
53
54 GridFunction nodes(&vfes);
55 mesh.SetNodalGridFunction(&nodes);
56 {
57 Array<int> dofs, vdofs;
58 GridFunction rdm(&vfes);
59 Vector h0(vfes.GetNDofs());
60 rdm.Randomize(seed);
61 rdm -= 0.5;
62 h0 = infinity();
63 for (int i = 0; i < mesh.GetNE(); i++)
64 {
65 vfes.GetElementDofs(i, dofs);
66 const double hi = mesh.GetElementSize(i);
67 for (int j = 0; j < dofs.Size(); j++)
68 {
69 h0(dofs[j]) = std::min(h0(dofs[j]), hi);
70 }
71 }
72 rdm.HostReadWrite();
73 for (int i = 0; i < vfes.GetNDofs(); i++)
74 {
75 for (int d = 0; d < dim; d++)
76 {
77 rdm(vfes.DofToVDof(i,d)) *= (0.25/p)*h0(i);
78 }
79 }
80 for (int i = 0; i < vfes.GetNBE(); i++)
81 {
82 vfes.GetBdrElementVDofs(i, vdofs);
83 for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
84 }
85 nodes -= rdm;
86 }
87
88 const Geometry::Type GeomType = mesh.GetElementBaseGeometry(0);
89 const IntegrationRule &ir = IntRules.Get(GeomType, 2*qpts-1);
90 const QuadratureInterpolator *sqi(sfes.GetQuadratureInterpolator(ir));
91 const QuadratureInterpolator *vqi(vfes.GetQuadratureInterpolator(ir));
92
93 const int NE(mesh.GetNE());
94 const int NQ(ir.GetNPoints());
95 const int ND(sfes.GetFE(0)->GetDof());
96 REQUIRE(ND == vfes.GetFE(0)->GetDof());
97
98 const ElementDofOrdering nat_ordering = ElementDofOrdering::NATIVE;
99 const ElementDofOrdering lex_ordering = ElementDofOrdering::LEXICOGRAPHIC;
100 const Operator *SRN(sfes.GetElementRestriction(nat_ordering));
101 const Operator *SRL(sfes.GetElementRestriction(lex_ordering));
102 const Operator *VRN(vfes.GetElementRestriction(nat_ordering));
103 const Operator *VRL(vfes.GetElementRestriction(lex_ordering));
104 MFEM_VERIFY(SRN, "No element sn-restriction operator found!");
105 MFEM_VERIFY(SRL, "No element sl-restriction operator found!");
106 MFEM_VERIFY(VRN, "No element vn-restriction operator found!");
107 MFEM_VERIFY(VRL, "No element vl-restriction operator found!");
108
109 const double rel_tol = 1e-12;
110
111 {
112 // Scalar
113 sqi->SetOutputLayout(q_layout);
114 Vector xe(1*ND*NE);
115 REQUIRE(xe.Size() == SRN->Height());
116 REQUIRE(SRN->Height() == SRL->Height());
117 // Full results
118 Vector sq_val_f(NQ*NE), sq_der_f(dim*NQ*NE), sq_pdr_f(dim*NQ*NE);
119 // Tensor results
120 Vector sq_val_t(NQ*NE), sq_der_t(dim*NQ*NE), sq_pdr_t(dim*NQ*NE);
121 {
122 // Full
123 SRN->Mult(x, xe);
124 sqi->DisableTensorProducts();
125
126 sqi->Values(xe, sq_val_f);
127
128 sqi->Derivatives(xe, sq_der_f);
129
130 sqi->PhysDerivatives(xe, sq_pdr_f);
131 }
132 {
133 // Tensor
134 SRL->Mult(x, xe);
135 sqi->EnableTensorProducts();
136
137 sqi->Values(xe, sq_val_t);
138
139 sqi->Derivatives(xe, sq_der_t);
140
141 sqi->PhysDerivatives(xe, sq_pdr_t);
142 }
143 double norm, rel_error;
144
145 norm = sq_val_f.Normlinf();
146 sq_val_f -= sq_val_t;
147 rel_error = sq_val_f.Normlinf()/norm;
148 if (verbose)
149 { std::cout << "sq_val rel. error = " << rel_error << std::endl; }
150 REQUIRE(rel_error <= rel_tol);
151
152 norm = sq_der_f.Normlinf();
153 sq_der_f -= sq_der_t;
154 rel_error = sq_der_f.Normlinf()/norm;
155 if (verbose)
156 { std::cout << "sq_der rel. error = " << rel_error << std::endl; }
157 REQUIRE(rel_error <= rel_tol);
158
159 norm = sq_pdr_f.Normlinf();
160 sq_pdr_f -= sq_pdr_t;
161 rel_error = sq_pdr_f.Normlinf()/norm;
162 if (verbose)
163 { std::cout << "sq_pdr rel. error = " << rel_error << std::endl; }
164 REQUIRE(rel_error <= rel_tol);
165 }
166
167 {
168 // Vector
169 vqi->SetOutputLayout(q_layout);
170 Vector ne(vdim*ND*NE);
171 REQUIRE(ne.Size() == VRN->Height());
172 REQUIRE(VRN->Height() == VRL->Height());
173 // Full results
174 Vector vq_val_f(dim*NQ*NE), vq_der_f(vdim*dim*NQ*NE),
175 vq_det_f(NQ*NE),
176 vq_pdr_f(vdim*dim*NQ*NE);
177 // Tensor results
178 Vector vq_val_t(dim*NQ*NE), vq_der_t(vdim*dim*NQ*NE),
179 vq_det_t(NQ*NE),
180 vq_pdr_t(vdim*dim*NQ*NE);
181 {
182 // Full
183 VRN->Mult(nodes, ne);
184 vqi->DisableTensorProducts();
185
186 vqi->Values(ne, vq_val_f);
187
188 vqi->Derivatives(ne, vq_der_f);
189
190 vqi->Determinants(ne, vq_det_f);
191
192 vqi->PhysDerivatives(ne, vq_pdr_f);
193 }
194 {
195 // Tensor
196 VRL->Mult(nodes, ne);
197 vqi->EnableTensorProducts();
198
199 vqi->Values(ne, vq_val_t);
200
201 vqi->Derivatives(ne, vq_der_t);
202
203 vqi->Determinants(ne, vq_det_t);
204
205 vqi->PhysDerivatives(ne, vq_pdr_t);
206 }
207 double norm, rel_error;
208
209 norm = vq_val_f.Normlinf();
210 vq_val_f -= vq_val_t;
211 rel_error = vq_val_f.Normlinf()/norm;
212 if (verbose)
213 { std::cout << "vq_val rel. error = " << rel_error << std::endl; }
214 REQUIRE(rel_error <= rel_tol);
215
216 norm = vq_der_f.Normlinf();
217 vq_der_f -= vq_der_t;
218 rel_error = vq_der_f.Normlinf()/norm;
219 if (verbose)
220 { std::cout << "vq_der rel. error = " << rel_error << std::endl; }
221 REQUIRE(rel_error <= rel_tol);
222
223 norm = vq_det_f.Normlinf();
224 vq_det_f -= vq_det_t;
225 rel_error = vq_det_f.Normlinf()/norm;
226 if (verbose)
227 { std::cout << "vq_det rel. error = " << rel_error << std::endl; }
228 REQUIRE(rel_error <= rel_tol);
229
230 norm = vq_pdr_f.Normlinf();
231 vq_pdr_f -= vq_pdr_t;
232 rel_error = vq_pdr_f.Normlinf()/norm;
233 if (verbose)
234 { std::cout << "vq_pdr rel. error = " << rel_error << std::endl; }
235 REQUIRE(rel_error <= rel_tol);
236 }
237
238 return true;
239 }
240
241 TEST_CASE("QuadratureInterpolator",
242 "[QuadratureInterpolator]"
243 "[CUDA]")
244 {
245 const auto d = GENERATE(2,3); // dimension
246 const auto p = GENERATE(range(1,7)); // element order, 1 <= p < 7
247 const auto q = GENERATE_COPY(p+1,p+2); // 1D quadrature points
248 const auto l = GENERATE(QVectorLayout::byNODES, QVectorLayout::byVDIM);
249 const auto nx = 3; // number of element in x
250 const auto ny = 3; // number of element in y
251 const auto nz = 3; // number of element in z
252 testQuadratureInterpolator(d, p, q, l, nx, ny, nz);
253 } // TEST_CASE "QuadratureInterpolator"
254