Difference between revisions of "Parametric domains"

From Medusa: Coordinate Free Mehless Method implementation
Jump to: navigation, search
Line 9: Line 9:
 
\end{align*}
 
\end{align*}
 
</math>
 
</math>
−
where $u(x,y)$ is the solution to the problem. Let's define $\Omega$ to be the interior of the parametrically given curve $f(t)$:
+
where $u(x,y)$ is the solution to the problem. Let's define $\Omega$ to be the interior of the parametrically given curve $r(t): [- 2 \pi, 0] \to \mathbb{R}^2$:
  
 
<math>
 
<math>
 
\begin{align*}
 
\begin{align*}
−
     r(t) &=& |\cos(1.5 t)| ^ {\sin(3t)}
+
     f(t) &= |\cos(-1.5 t)| ^ {\sin(-3t)} \\
−
         f(t) &=& (r \cos(t), r \sin(t))
+
         r(t) &= (r \cos(-t), r \sin(-t))
 
\end{align*}
 
\end{align*}
 
</math>
 
</math>
 +
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:
 +
<syntaxhighlight lang="c++" line>
 +
auto example_r = [](Vec1d t) {
 +
    t(0) = -t(0);
 +
    double f = pow(abs(cos(1.5 * t(0))), sin(3 * t(0)));
 +
    return Vec2d(f * cos(t(0)), f * sin(t(0)));
 +
};
  
−
See [[Positioning of computational nodes]] TODO.
+
auto der_example_r = [](Vec1d t) {
 +
    t(0) = -t(0);
 +
    double f = pow(abs(cos(1.5 * t(0))), sin(3 * t(0)));
 +
    double der_f = (-1.5 * pow(abs(cos(1.5 * t(0))),
 +
            sin(3 * t(0))) * sin(3 * t(0)) * sin(1.5 * t(0))  +  3 * pow(abs(cos(1.5 * t(0))),
 +
                    sin(3 * t(0))) * cos(3 * t(0)) * cos(1.5 * t(0)) *
 +
                            log(abs(cos(1.5 * t(0))))) / cos(1.5 * t(0));
 +
 
 +
    Eigen::Matrix<double, 2, 1> jm;
 +
    jm.col(0) << -(der_f * cos(t(0)) - f * sin(t(0))), -(der_f * sin(t(0)) + f * cos(t(0)));
 +
 
 +
    return jm;
 +
};
 +
 
 +
// Define parametric curve's domain.
 +
BoxShape<Vec1d> bs(Vec<double, 1>{- 2 * PI}, Vec<double, 1>{0.0});
 +
DomainDiscretization<Vec1d> param_d(bs);
 +
</syntaxhighlight>
 +
 
 +
We implemented $r$ and $\textbf{J}$ as lambda expressions and $r$'s domain as an instance of <code>DomainDiscretization</code> class. Now we can fill the target domain with nodes given by $r$ and the density function $h$ by using <code>GeneralSurfaceFill</code>. For the sake of example, we choose a gradient-like density function.
 +
<syntaxhighlight lang="c++" line>
 +
UnknownShape<Vec2d> shape;
 +
DomainDiscretization<Vec2d> domain(shape);
 +
 
 +
auto gradient_h = [](Vec2d p){
 +
    double h_0 = 0.005;
 +
    double h_m = 0.03 - h_0;
 +
 
 +
    return  (0.5 * h_m * (p(0) + p(1) + 3.0) + h_0) / 5.0;
 +
};
 +
 
 +
GeneralSurfaceFill<Vec2d, Vec1d> gsf;
 +
domain.fill(gsf, param_d, example_r, der_example_r, gradient_h);
 +
</syntaxhighlight>
 +
 
 +
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.
 +
<syntaxhighlight lang="c++" line>
 +
GeneralFill<Vec2d> gf;
 +
domain.fill(gf, gradient_h);
 +
 
 +
int N = domain.size();
 +
domain.findSupport(FindClosest(9));
 +
</syntaxhighlight>
 +
 
 +
Finally, we translate the equations of our problem into code and solve the problem, as in other examples.
 +
<syntaxhighlight lang="c++" line>
 +
int m = 2;  // basis order
 +
Monomials<Vec2d> mon(m);
 +
RBFFD<Polyharmonic<double, 3>, Vec2d> approx({}, mon);
 +
 
 +
// Compute the shapes (we only need the Laplacian).
 +
auto storage = domain.computeShapes<sh::lap>(approx);
 +
 
 +
Eigen::SparseMatrix<double, Eigen::RowMajor> M(N, N);
 +
Eigen::VectorXd rhs(N); rhs.setZero();
 +
M.reserve(storage.supportSizes());
 +
 
 +
// Construct implicit operators over our storage.
 +
auto op = storage.implicitOperators(M, rhs);
 +
 
 +
for (int i : domain.interior()) {
 +
    op.lap(i) = 0.5; // set the case for nodes in the domain
 +
}
 +
for (int i : domain.boundary()) {
 +
    op.value(i) = 0.0; // enforce the boundary conditions
 +
}
 +
 
 +
Eigen::BiCGSTAB<decltype(M), Eigen::IncompleteLUT<double>> solver;
 +
solver.compute(M);
 +
ScalarFieldd u = solver.solve(rhs);
 +
</syntaxhighlight>
 +
 
 +
Here is the plot of $u(x, y)$.
 +
 
 +
TODO
 +
 
 +
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].
 +
 
 +
See [[Positioning of computational nodes]] (TODO) for the theory behind the node placing algorithm.

Revision as of 15:49, 8 August 2019

Go back to Examples.

Variable node density and dirchlet 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*} \] where $u(x,y)$ is the solution to the problem. Let's define $\Omega$ to be the interior of the parametrically given curve $r(t): [- 2 \pi, 0] \to \mathbb{R}^2$\[ \begin{align*} f(t) &= |\cos(-1.5 t)| ^ {\sin(-3t)} \\ r(t) &= (r \cos(-t), r \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 = [](Vec1d t) {
 2     t(0) = -t(0);
 3     double f = pow(abs(cos(1.5 * t(0))), sin(3 * t(0)));
 4     return Vec2d(f * cos(t(0)), f * sin(t(0)));
 5 };
 6 
 7 auto der_example_r = [](Vec1d t) {
 8     t(0) = -t(0);
 9     double f = pow(abs(cos(1.5 * t(0))), sin(3 * t(0)));
10     double der_f = (-1.5 * pow(abs(cos(1.5 * t(0))),
11             sin(3 * t(0))) * sin(3 * t(0)) * sin(1.5 * t(0))  +  3 * pow(abs(cos(1.5 * t(0))),
12                     sin(3 * t(0))) * cos(3 * t(0)) * cos(1.5 * t(0)) *
13                             log(abs(cos(1.5 * t(0))))) / cos(1.5 * t(0));
14 
15     Eigen::Matrix<double, 2, 1> jm;
16     jm.col(0) << -(der_f * cos(t(0)) - f * sin(t(0))), -(der_f * sin(t(0)) + f * cos(t(0)));
17 
18     return jm;
19 };
20 
21 // Define parametric curve's domain.
22 BoxShape<Vec1d> bs(Vec<double, 1>{- 2 * PI}, Vec<double, 1>{0.0});
23 DomainDiscretization<Vec1d> param_d(bs);

We implemented $r$ and $\textbf{J}$ as lambda expressions and $r$'s domain as an instance of DomainDiscretization 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 example, we choose a gradient-like density function. 

 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_d, 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 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)$.

TODO

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.

See Positioning of computational nodes (TODO) for the theory behind the node placing algorithm.

Retrieved from "http://e6.ijs.si/medusa/wiki/index.php?title=Parametric_domains&oldid=2685"