MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
sbm_aux.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 
14 using namespace std;
15 using namespace mfem;
16 
17 double point_inside_trigon(const Vector px, Vector p1, Vector p2, Vector p3)
18 {
19  Vector v0 = p1;
20  Vector v1 = p2; v1 -=p1;
21  Vector v2 = p3; v2 -=p1;
22  double p, q;
23  p = ((px(0)*v2(1)-px(1)*v2(0))-(v0(0)*v2(1)-v0(1)*v2(0))) /
24  (v1(0)*v2(1)-v1(1)*v2(0));
25  q = -((px(0)*v1(1)-px(1)*v1(0))-(v0(0)*v1(1)-v0(1)*v1(0))) /
26  (v1(0)*v2(1)-v1(1)*v2(0));
27 
28  return (p > 0 && q > 0 && 1-p-q > 0) ? -1.0 : 1.0;
29 }
30 
31 // 1 is inside the doughnut, -1 is outside.
32 double doughnut_cheese(const Vector &coord)
33 {
34  // map [0,1] to [-1,1].
35  double x = 2*coord(0)-1.0, y = 2*coord(1)-1.0, z = 2*coord(2)-1.0;
36 
37  bool doughnut;
38  const double R = 0.8, r = 0.15;
39  const double t = R - std::sqrt(x*x + y*y);
40  doughnut = t*t + z*z - r*r <= 0;
41 
42  bool cheese;
43  x = 3.0*x, y = 3.0*y, z = 3.0*z;
44  cheese = (x*x + y*y - 4.0) * (x*x + y*y - 4.0) +
45  (z*z - 1.0) * (z*z - 1.0) +
46  (y*y + z*z - 4.0) * (y*y + z*z - 4.0) +
47  (x*x - 1.0) * (x*x - 1.0) +
48  (z*z + x*x - 4.0) * (z*z + x*x - 4.0) +
49  (y*y - 1.0) * (y*y - 1.0) - 15.0 <= 0.0;
50 
51  return (doughnut || cheese) ? 1.0 : -1.0;
52 }
53 
54 /// Analytic distance to the 0 level set. Positive value if the point is inside
55 /// the domain, and negative value if outside.
56 double dist_value(const Vector &x, const int type)
57 {
58  if (type == 1 || type == 2) // circle of radius 0.2 - centered at 0.5, 0.5
59  {
60  const double ring_radius = 0.2;
61  Vector xc(x.Size());
62  xc = 0.5;
63  xc -= x;
64  return xc.Norml2() - ring_radius; // positive is the domain
65  }
66  else if (type == 3) // walls at y = 0.0
67  {
68  return x(1);
69  }
70  else if (type == 4)
71  {
72  const int num_circ = 3;
73  double rad[num_circ] = {0.3, 0.15, 0.2};
74  double c[num_circ][2] = { {0.6, 0.6}, {0.3, 0.3}, {0.25, 0.75} };
75 
76  const double xc = x(0), yc = x(1);
77 
78  // circle 0
79  double r0 = (xc-c[0][0])*(xc-c[0][0]) + (yc-c[0][1])*(yc-c[0][1]);
80  r0 = (r0 > 0) ? std::sqrt(r0) : 0.0;
81  if (r0 <= 0.2) { return -1.0; }
82 
83  for (int i = 0; i < num_circ; i++)
84  {
85  double r = (xc-c[i][0])*(xc-c[i][0]) + (yc-c[i][1])*(yc-c[i][1]);
86  r = (r > 0) ? std::sqrt(r) : 0.0;
87  if (r <= rad[i]) { return 1.0; }
88  }
89 
90  // rectangle 1
91  if (0.7 <= xc && xc <= 0.8 && 0.1 <= yc && yc <= 0.8) { return 1.0; }
92 
93  // rectangle 2
94  if (0.3 <= xc && xc <= 0.8 && 0.15 <= yc && yc <= 0.2) { return 1.0; }
95  return -1.0;
96  }
97  else if (type == 5) // square of side 0.2 centered at 0.75, 0.25
98  {
99  double square_side = 0.2;
100  Vector xc(x.Size());
101  xc = 0.75; xc(1) = 0.25;
102  xc -= x;
103  if (abs(xc(0)) > 0.5*square_side || abs(xc(1)) > 0.5*square_side)
104  {
105  return 1.0;
106  }
107  else
108  {
109  return -1.0;
110  }
111  return 0.0;
112  }
113  else if (type == 6) // Triangle
114  {
115  Vector p1(x.Size()), p2(x.Size()), p3(x.Size());
116  p1(0) = 0.25; p1(1) = 0.4;
117  p2(0) = 0.1; p2(1) = 0.1;
118  p3(0) = 0.4; p3(1) = 0.1;
119  return point_inside_trigon(x, p1, p2, p3);
120  }
121  else if (type == 7) // circle of radius 0.2 - centered at 0.5, 0.6
122  {
123  Vector xc(x.Size());
124  xc = 0.5; xc(1) = 0.6;
125  xc -= x;
126  return xc.Norml2() - 0.2;
127  }
128  else if (type == 8) { return doughnut_cheese(x); }
129  else
130  {
131  MFEM_ABORT(" Function type not implement yet.");
132  }
133  return 0.;
134 }
135 
136 /// Level set coefficient: +1 inside the true domain, -1 outside.
138 {
139 private:
140  int type;
141 
142 public:
144  : Coefficient(), type(type_) { }
145 
146  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
147  {
148  Vector x(3);
149  T.Transform(ip, x);
150  double dist = dist_value(x, type);
151  return (dist >= 0.0) ? 1.0 : -1.0;
152  }
153 };
154 
155 /// Combination of level sets: +1 inside the true domain, -1 outside.
157 {
158 private:
160 
161 public:
163 
165  { dls.Append(&dls_); }
166 
167  int GetNLevelSets() { return dls.Size(); }
168 
169  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
170  {
171  MFEM_VERIFY(dls.Size() > 0,
172  "Add at least 1 Dist_level_Set_Coefficient to the Combo.");
173  double dist = dls[0]->Eval(T, ip);
174  for (int j = 1; j < dls.Size(); j++)
175  {
176  dist = min(dist, dls[j]->Eval(T, ip));
177  }
178  return (dist >= 0.0) ? 1.0 : -1.0;
179  }
180 };
181 
182 /// Distance vector to the zero level-set.
184 {
185 private:
186  int type;
187 
188 public:
189  Dist_Vector_Coefficient(int dim_, int type_)
190  : VectorCoefficient(dim_), type(type_) { }
191 
192  using VectorCoefficient::Eval;
193 
194  virtual void Eval(Vector &p, ElementTransformation &T,
195  const IntegrationPoint &ip)
196  {
197  Vector x;
198  T.Transform(ip, x);
199  const int dim = x.Size();
200  p.SetSize(dim);
201  if (type == 1 || type == 2)
202  {
203  double dist0 = dist_value(x, type);
204  for (int i = 0; i < dim; i++) { p(i) = 0.5 - x(i); }
205  double length = p.Norml2();
206  p *= dist0/length;
207  }
208  else if (type == 3)
209  {
210  double dist0 = dist_value(x, type);
211  p(0) = 0.;
212  p(1) = -dist0;
213  }
214  }
215 };
216 
217 /// Boundary conditions - Dirichlet
218 double homogeneous(const Vector &x)
219 {
220  return 0.0;
221 }
222 
224 {
225  double xy_p = 2.; // exponent for level set 2 where u = x^p+y^p;
226  return pow(x(0), xy_p) + pow(x(1), xy_p);
227 }
228 
230 {
231  return 1./(M_PI*M_PI)*std::sin(M_PI*x(0)*x(1));
232 }
233 
234 /// Boundary conditions - Neumann
235 /// Normal vector for level_set_type = 1. Circle centered at [0.5 , 0.5]
236 void normal_vector_1(const Vector &x, Vector &p)
237 {
238  p.SetSize(x.Size());
239  p(0) = x(0)-0.5;
240  p(1) = x(1)-0.5; // center of circle at [0.5, 0.5]
241  p /= p.Norml2();
242  p *= -1;
243 }
244 
245 /// Normal vector for level_set_type = 7. Circle centered at [0.5 , 0.6]
246 void normal_vector_2(const Vector &x, Vector &p)
247 {
248  p.SetSize(x.Size());
249  p(0) = x(0)-0.5;
250  p(1) = x(1)-0.6; // center of circle at [0.5, 0.6]
251  p /= p.Norml2();
252  p *= -1;
253 }
254 
255 /// Neumann condition for exponent based solution
256 double traction_xy_exponent(const Vector &x)
257 {
258  double xy_p = 2;
259  Vector gradient(2);
260  gradient(0) = xy_p*x(0);
261  gradient(1) = xy_p*x(1);
262  Vector normal(2);
263  normal_vector_1(x, normal);
264  return 1.0*(gradient*normal);
265 }
266 
267 /// `f` for the Poisson problem (-nabla^2 u = f).
268 double rhs_fun_circle(const Vector &x)
269 {
270  return 1;
271 }
272 
273 double rhs_fun_xy_exponent(const Vector &x)
274 {
275  double xy_p = 2.; // exponent for level set 2 where u = x^p+y^p;
276  double coeff = std::max(xy_p*(xy_p-1), 1.);
277  double expon = std::max(0., xy_p-2);
278  if (xy_p == 1)
279  {
280  return 0.;
281  }
282  else
283  {
284  return -coeff*std::pow(x(0), expon) - coeff*std::pow(x(1), expon);
285  }
286 }
287 
289 {
290  return std::sin(M_PI*x(0)*x(1))*(x(0)*x(0)+x(1)*x(1));
291 }
double rhs_fun_circle(const Vector &x)
f for the Poisson problem (-nabla^2 u = f).
Definition: sbm_aux.hpp:268
Base class for vector Coefficients that optionally depend on time and space.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:812
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
double dirichlet_velocity_xy_exponent(const Vector &x)
Definition: sbm_aux.hpp:223
double homogeneous(const Vector &x)
Boundary conditions - Dirichlet.
Definition: sbm_aux.hpp:218
virtual void Eval(Vector &p, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Definition: sbm_aux.hpp:194
Distance vector to the zero level-set.
Definition: sbm_aux.hpp:183
Level set coefficient: +1 inside the true domain, -1 outside.
Definition: sbm_aux.hpp:137
Dist_Vector_Coefficient(int dim_, int type_)
Definition: sbm_aux.hpp:189
double point_inside_trigon(const Vector px, Vector p1, Vector p2, Vector p3)
Definition: sbm_aux.hpp:17
double rhs_fun_xy_sinusoidal(const Vector &x)
Definition: sbm_aux.hpp:288
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Definition: sbm_aux.hpp:169
Combination of level sets: +1 inside the true domain, -1 outside.
Definition: sbm_aux.hpp:156
void Add_Level_Set_Coefficient(Dist_Level_Set_Coefficient &dls_)
Definition: sbm_aux.hpp:164
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Definition: sbm_aux.hpp:146
Dist_Level_Set_Coefficient(int type_)
Definition: sbm_aux.hpp:143
double rhs_fun_xy_exponent(const Vector &x)
Definition: sbm_aux.hpp:273
Class for integration point with weight.
Definition: intrules.hpp:25
double dirichlet_velocity_xy_sinusoidal(const Vector &x)
Definition: sbm_aux.hpp:229
void normal_vector_1(const Vector &x, Vector &p)
Definition: sbm_aux.hpp:236
int dim
Definition: ex24.cpp:53
void normal_vector_2(const Vector &x, Vector &p)
Normal vector for level_set_type = 7. Circle centered at [0.5 , 0.6].
Definition: sbm_aux.hpp:246
RefCoord t[3]
double traction_xy_exponent(const Vector &x)
Neumann condition for exponent based solution.
Definition: sbm_aux.hpp:256
Vector data type.
Definition: vector.hpp:60
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
double doughnut_cheese(const Vector &coord)
Definition: sbm_aux.hpp:32
double dist_value(const Vector &x, const int type)
Definition: sbm_aux.hpp:56