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 "catch.hpp"
14
15 using namespace mfem;
16
17 namespace surf_blf
18 {
19
20 Mesh * GetCylMesh(int type);
21
22 void trans_skew_cyl(const Vector &x, Vector &r);
23
24 void sigmaFunc(const Vector &x, DenseMatrix &s);
25
uExact(const Vector & x)26 double uExact(const Vector &x)
27 {
28 return (0.25 * (2.0 + x[0]) - x[2]) * (x[2] + 0.25 * (2.0 + x[0]));
29 }
30
31 TEST_CASE("Embedded Surface Diffusion",
32 "[DiffusionIntegrator]")
33 {
34 for (int type = (int) Element::TRIANGLE;
35 type <= (int) Element::QUADRILATERAL;
36 type++)
37 {
38 int order = 3;
39
40 Mesh *mesh = GetCylMesh(type);
41 int dim = mesh->Dimension();
42
43 mesh->SetCurvature(3);
44 mesh->Transform(trans_skew_cyl);
45
46 H1_FECollection fec(order, dim);
47 FiniteElementSpace fespace(mesh, &fec);
48
49 Array<int> ess_tdof_list;
50 if (mesh->bdr_attributes.Size())
51 {
52 Array<int> ess_bdr(mesh->bdr_attributes.Max());
53 ess_bdr = 1;
54 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
55 }
56
57 LinearForm b(&fespace);
58 ConstantCoefficient one(1.0);
59 b.AddDomainIntegrator(new DomainLFIntegrator(one));
60 b.Assemble();
61
62 GridFunction x(&fespace);
63 x = 0.0;
64
65 BilinearForm a(&fespace);
66 MatrixFunctionCoefficient sigma(3, sigmaFunc);
67 a.AddDomainIntegrator(new DiffusionIntegrator(sigma));
68 a.Assemble();
69
70 OperatorPtr A;
71 Vector B, X;
72 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
73
74 GSSmoother M((SparseMatrix&)(*A));
75 PCG(*A, M, B, X, 0, 200, 1e-12, 0.0);
76
77 a.RecoverFEMSolution(X, b, x);
78
79 FunctionCoefficient uCoef(uExact);
80 double err = x.ComputeL2Error(uCoef);
81
82 REQUIRE(err < 1.5e-3);
83
84 delete mesh;
85 }
86 }
87
GetCylMesh(int type)88 Mesh * GetCylMesh(int type)
89 {
90 Mesh * mesh = NULL;
91
92 if (type == (int) Element::TRIANGLE)
93 {
94 mesh = new Mesh(2, 12, 16, 8, 3);
95
96 mesh->AddVertex(-1.0, -1.0, 0.0);
97 mesh->AddVertex( 1.0, -1.0, 0.0);
98 mesh->AddVertex( 1.0, 1.0, 0.0);
99 mesh->AddVertex(-1.0, 1.0, 0.0);
100 mesh->AddVertex(-1.0, -1.0, 1.0);
101 mesh->AddVertex( 1.0, -1.0, 1.0);
102 mesh->AddVertex( 1.0, 1.0, 1.0);
103 mesh->AddVertex(-1.0, 1.0, 1.0);
104 mesh->AddVertex( 0.0, -1.0, 0.5);
105 mesh->AddVertex( 1.0, 0.0, 0.5);
106 mesh->AddVertex( 0.0, 1.0, 0.5);
107 mesh->AddVertex(-1.0, 0.0, 0.5);
108
109 mesh->AddTriangle(0, 1, 8);
110 mesh->AddTriangle(1, 5, 8);
111 mesh->AddTriangle(5, 4, 8);
112 mesh->AddTriangle(4, 0, 8);
113 mesh->AddTriangle(1, 2, 9);
114 mesh->AddTriangle(2, 6, 9);
115 mesh->AddTriangle(6, 5, 9);
116 mesh->AddTriangle(5, 1, 9);
117 mesh->AddTriangle(2, 3, 10);
118 mesh->AddTriangle(3, 7, 10);
119 mesh->AddTriangle(7, 6, 10);
120 mesh->AddTriangle(6, 2, 10);
121 mesh->AddTriangle(3, 0, 11);
122 mesh->AddTriangle(0, 4, 11);
123 mesh->AddTriangle(4, 7, 11);
124 mesh->AddTriangle(7, 3, 11);
125
126 mesh->AddBdrSegment(0, 1, 1);
127 mesh->AddBdrSegment(1, 2, 1);
128 mesh->AddBdrSegment(2, 3, 1);
129 mesh->AddBdrSegment(3, 0, 1);
130 mesh->AddBdrSegment(5, 4, 2);
131 mesh->AddBdrSegment(6, 5, 2);
132 mesh->AddBdrSegment(7, 6, 2);
133 mesh->AddBdrSegment(4, 7, 2);
134 }
135 else if (type == (int) Element::QUADRILATERAL)
136 {
137 mesh = new Mesh(2, 8, 4, 8, 3);
138
139 mesh->AddVertex(-1.0, -1.0, 0.0);
140 mesh->AddVertex( 1.0, -1.0, 0.0);
141 mesh->AddVertex( 1.0, 1.0, 0.0);
142 mesh->AddVertex(-1.0, 1.0, 0.0);
143 mesh->AddVertex(-1.0, -1.0, 1.0);
144 mesh->AddVertex( 1.0, -1.0, 1.0);
145 mesh->AddVertex( 1.0, 1.0, 1.0);
146 mesh->AddVertex(-1.0, 1.0, 1.0);
147
148 mesh->AddQuad(0, 1, 5, 4);
149 mesh->AddQuad(1, 2, 6, 5);
150 mesh->AddQuad(2, 3, 7, 6);
151 mesh->AddQuad(3, 0, 4, 7);
152
153 mesh->AddBdrSegment(0, 1, 1);
154 mesh->AddBdrSegment(1, 2, 1);
155 mesh->AddBdrSegment(2, 3, 1);
156 mesh->AddBdrSegment(3, 0, 1);
157 mesh->AddBdrSegment(5, 4, 2);
158 mesh->AddBdrSegment(6, 5, 2);
159 mesh->AddBdrSegment(7, 6, 2);
160 mesh->AddBdrSegment(4, 7, 2);
161 }
162 else
163 {
164 MFEM_ABORT("Unrecognized mesh type " << type << "!");
165 }
166 mesh->FinalizeTopology();
167
168 return mesh;
169 }
170
trans_skew_cyl(const Vector & x,Vector & r)171 void trans_skew_cyl(const Vector &x, Vector &r)
172 {
173 r.SetSize(3);
174
175 double tol = 1e-6;
176 double theta = 0.0;
177 if (fabs(x[1] + 1.0) < tol)
178 {
179 theta = 0.25 * M_PI * (x[0] - 2.0);
180 }
181 else if (fabs(x[0] - 1.0) < tol)
182 {
183 theta = 0.25 * M_PI * x[1];
184 }
185 else if (fabs(x[1] - 1.0) < tol)
186 {
187 theta = 0.25 * M_PI * (2.0 - x[0]);
188 }
189 else if (fabs(x[0] + 1.0) < tol)
190 {
191 theta = 0.25 * M_PI * (4.0 - x[1]);
192 }
193 else
194 {
195 mfem::out << "side not recognized "
196 << x[0] << " " << x[1] << " " << x[2] << std::endl;
197 }
198
199 r[0] = cos(theta);
200 r[1] = sin(theta);
201 r[2] = 0.25 * (2.0 * x[2] - 1.0) * (r[0] + 2.0);
202 }
203
sigmaFunc(const Vector & x,DenseMatrix & s)204 void sigmaFunc(const Vector &x, DenseMatrix &s)
205 {
206 s.SetSize(3);
207 double a = 17.0 - 2.0 * x[0] * (1.0 + x[0]);
208 s(0,0) = 0.5 + x[0] * x[0] * (8.0 / a - 0.5);
209 s(0,1) = x[0] * x[1] * (8.0 / a - 0.5);
210 s(0,2) = 0.0;
211 s(1,0) = s(0,1);
212 s(1,1) = 0.5 * x[0] * x[0] + 8.0 * x[1] * x[1] / a;
213 s(1,2) = 0.0;
214 s(2,0) = 0.0;
215 s(2,1) = 0.0;
216 s(2,2) = a / 32.0;
217 }
218
219 } // namespace surf_blf
220