Difference between revisions of "Parametric domains"
| (One intermediate revision by the same user not shown) | |||
| Line 103: | Line 103: | ||
[[File:Poisson_2D.png|800px]] | [[File:Poisson_2D.png|800px]] | ||
| − | The whole example can be found as [https://gitlab.com/e62Lab/medusa parametric_domain_2D.cpp] along with the plot script that can be run by Matlab or Octave [https://gitlab.com/e62Lab/medusa parametric_domain_2D.m | + | The whole example can be found as [https://gitlab.com/e62Lab/medusa/-/blob/dev/examples/poisson_equation/parametric_domain_2D.cpp parametric_domain_2D.cpp] along with the plot script that can be run by Matlab or Octave [https://gitlab.com/e62Lab/medusa/-/blob/dev/examples/poisson_equation/parametric_domain_2D.m parametric_domain_2D.m]. |
| − | |||
| − | |||
Latest revision as of 12:26, 1 September 2020
Go back to Examples.
Medusa supports creating domains bounded by parametric curves and (hyper-)surfaces and filling them and their interiors with node distributions of arbitrary densities.
See Positioning of computational nodes for details on node positioning algorithms.
Variable node density and Dirichlet boundary conditions in 2D
With Medusa, we can also solve partial differential equations on parametric domains. Consider the solution of a simple 2D Poisson equation with Dirichlet boundary conditions\[ \begin{align*} \Delta u &= 0.5 &&\text{in } \Omega, \\ u &= 0 &&\text{on } \partial \Omega, \end{align*} \] Let's define $\Omega$ to be the interior of the parametrically given curve $r(t): [0, 2 \pi] \to \mathbb{R}^2$\[ \begin{align*} f(t) &= |\cos(1.5 t)| ^ {\sin(3t)} \\ r(t) &= (f(t) \cos(t), f(t) \sin(t)) \end{align*} \] Using tools such as Wolfram Mathematica, we can find the Jacobian matrix ($\textbf{J}$) of $r$. Now, we can translate the definition of $r$ into code:
1 auto example_r = [](Vec<double, 1> t) {
2 double r = pow(abs(cos(1.5 * t(0))), sin(3 * t(0)));
3 return Vec2d(r * cos(t(0)), r * sin(t(0)));
4 };
5
6 auto der_example_r = [](Vec<double, 1> t) {
7 double r = pow(abs(cos(1.5 * t(0))), sin(3 * t(0)));
8 double der_r = (-1.5 * pow(abs(cos(1.5 * t(0))),
9 sin(3 * t(0))) * sin(3 * t(0)) * sin(1.5 * t(0)) +
10 3 * pow(abs(cos(1.5 * t(0))),
11 sin(3 * t(0))) * cos(3 * t(0)) * cos(1.5 * t(0))
12 * log(abs(cos(1.5 * t(0))))) / cos(1.5 * t(0));
13
14 Eigen::Matrix<double, 2, 1> jm;
15 jm.col(0) << der_r * cos(t(0)) - r * sin(t(0)), der_r * sin(t(0)) + r * cos(t(0));
16
17 return jm;
18 };
19
20 // Define parametric curve's domain.
21 BoxShape<Vec1d> param_bs(Vec<double, 1>{0.0}, Vec<double, 1>{2 * PI});
We implemented $r$ and $\textbf{J}$ as lambda expressions and $r$'s domain as an instance of BoxShape class. Now we can fill the target domain with nodes given by $r$ and the density function $h$ by using GeneralSurfaceFill. For the sake of this example, we choose a gradient-like density function gradient_h.
1 UnknownShape<Vec2d> shape;
2 DomainDiscretization<Vec2d> domain(shape);
3
4 auto gradient_h = [](Vec2d p){
5 double h_0 = 0.005;
6 double h_m = 0.03 - h_0;
7
8 return (0.5 * h_m * (p(0) + p(1) + 3.0) + h_0) / 5.0;
9 };
10
11 GeneralSurfaceFill<Vec2d, Vec1d> gsf;
12 domain.fill(gsf, param_bs, example_r, der_example_r, gradient_h);
After that, we can fill the interior of the domain in the usual way and choose the closest 9 nodes of each node as its support.
1 GeneralFill<Vec2d> gf;
2 domain.fill(gf, gradient_h);
3
4 int N = domain.size();
5 domain.findSupport(FindClosest(9));
Finally, we translate the partial differential equations of our problem into code and solve the problem, as in other examples.
1 int m = 2; // basis order
2 Monomials<Vec2d> mon(m);
3 RBFFD<Polyharmonic<double, 3>, Vec2d> approx({}, mon);
4
5 // Compute the shapes (we only need the Laplacian).
6 auto storage = domain.computeShapes<sh::lap>(approx);
7
8 Eigen::SparseMatrix<double, Eigen::RowMajor> M(N, N);
9 Eigen::VectorXd rhs(N); rhs.setZero();
10 M.reserve(storage.supportSizes());
11
12 // Construct implicit operators over our storage.
13 auto op = storage.implicitOperators(M, rhs);
14
15 for (int i : domain.interior()) {
16 op.lap(i) = 0.5; // set the case for nodes in the domain
17 }
18 for (int i : domain.boundary()) {
19 op.value(i) = 0.0; // enforce the boundary conditions
20 }
21
22 Eigen::BiCGSTAB<decltype(M), Eigen::IncompleteLUT<double>> solver;
23 solver.compute(M);
24 ScalarFieldd u = solver.solve(rhs);
Here is the plot of $u(x, y)$.
The whole example can be found as parametric_domain_2D.cpp along with the plot script that can be run by Matlab or Octave parametric_domain_2D.m.
