Coupled domains - Revision history

Jureslak at 08:31, 18 May 2019

2019-05-18T08:31:25Z

Jureslak at 18:50, 18 April 2019

2019-04-18T18:50:22Z

Jureslak at 18:49, 18 April 2019

2019-04-18T18:49:52Z

Jureslak at 18:48, 18 April 2019

2019-04-18T18:48:38Z

Jureslak at 18:47, 18 April 2019

2019-04-18T18:47:41Z

Jureslak: /* Discretisation */

2019-04-18T18:46:36Z

‎Discretisation

Jureslak at 18:41, 18 April 2019

2019-04-18T18:41:13Z

Jureslak at 18:40, 18 April 2019

2019-04-18T18:40:57Z

Jureslak at 18:40, 18 April 2019

2019-04-18T18:40:01Z

BlazS at 17:02, 18 April 2019

2019-04-18T17:02:57Z

← Older revision		Revision as of 08:31, 18 May 2019
Line 1:		Line 1:
		+	Go back to [[Medusa#Examples\|Examples]].
		+
	In this example we will show how one can solve coupled problems with the Medusa library. We will solve the stationary heat equation on concentric circles where the thermal conductivity $\lambda$ is different for each circle. The inner circle will be heated up, while the outside temperature will be held at zero.		In this example we will show how one can solve coupled problems with the Medusa library. We will solve the stationary heat equation on concentric circles where the thermal conductivity $\lambda$ is different for each circle. The inner circle will be heated up, while the outside temperature will be held at zero.

Line 98:		Line 100:

	The example can be found here [https://gitlab.com/e62Lab/medusa/blob/dev/examples/coupled_domains/poisson_coupled_domains.cpp poisson_coupled_domains.cpp], and the matlab plot script as [https://gitlab.com/e62Lab/medusa/blob/dev/examples/coupled_domains/poisson_coupled_domains.m poisson_coupled_domains.m].		The example can be found here [https://gitlab.com/e62Lab/medusa/blob/dev/examples/coupled_domains/poisson_coupled_domains.cpp poisson_coupled_domains.cpp], and the matlab plot script as [https://gitlab.com/e62Lab/medusa/blob/dev/examples/coupled_domains/poisson_coupled_domains.m poisson_coupled_domains.m].
		+
		+	Go back to [[Medusa#Examples\|Examples]].

@@ Line 12: / Line 12: @@
      u = 0 \qquad &\text{on} \quad \partial \Omega_1 \\
      u = v \qquad &\text{on} \quad \partial \Omega_2 \\
-     \lambda_2 \dpar{u}{n} + \lambda_1 \dpar{v}{n} = 0 \qquad &\text{on} \quad \partial \Omega_2. \\
+     \lambda_2 \frac{\partial u}{\partial n} + \lambda_1 \frac{\partial v}{\partial n} = 0 \qquad &\text{on} \quad \partial \Omega_2. \\
 \end{align*}
 We will see that the task is quite simple using Medusa.

@@ Line 12: / Line 12: @@
      u = 0 \qquad &\text{on} \quad \partial \Omega_1 \\
      u = v \qquad &\text{on} \quad \partial \Omega_2 \\
-     \lambda_2 \frac{\partial u}{\partial n} + \lambda_1 \frac{\partial v}{\partial n} = 0 \qquad &\text{on} \quad \partial \Omega_2. \\
+     \lambda_2 \dpar{u}{n} + \lambda_1 \dpar{v}{n} = 0 \qquad &\text{on} \quad \partial \Omega_2. \\
 \end{align*}
 We will see that the task is quite simple using Medusa.

@@ Line 12: / Line 12: @@
      u = 0 \qquad &\text{on} \quad \partial \Omega_1 \\
      u = v \qquad &\text{on} \quad \partial \Omega_2 \\
-     \lambda_2 \frac{\partial u}{\partial n_1} + \lambda_1 \frac{\partial v}{\partial n_2} = 0 \qquad &\text{on} \quad \partial \Omega_2. \\
+     \lambda_2 \frac{\partial u}{\partial n} + \lambda_1 \frac{\partial v}{\partial n} = 0 \qquad &\text{on} \quad \partial \Omega_2. \\
 \end{align*}
 We will see that the task is quite simple using Medusa.

@@ Line 44: / Line 44: @@
      inner.fill(fill_engine, h); // this is the inner circle
 </syntaxhighlight>
-The snippet above shows how to add the the nodes on the common boundary to the inner domain. The normals of the shared boundary nodes belonging to the inner domain must point into the opposite direction as the normals belonging to the shared boundary nodes in the outer domain. Additionally, mapping of node indices from outer domain to the index of the same node in the inner domain was made. After this is done, the insides of both domains are filled. The size of the matrix $M$ of the system must be equal to the sum of numbers of nodes in the inner and outer domains.
+The snippet above shows how to add the the nodes on the common boundary to the inner domain. The normals of the shared boundary nodes belonging to the inner domain must point into the opposite direction as the normals belonging to the shared boundary nodes in the outer domain. Additionally, mapping of node indices from outer domain to the index of the same node in the inner domain was made. After this is done, the insides of both domains are filled.
 <syntaxhighlight lang="c++" line>
      Eigen::SparseMatrix<double, Eigen::RowMajor> M(N_inner+N_outer, N_inner+N_outer);

@@ Line 31: / Line 31: @@
      auto outer_bnd = outer.positions().filter([=](const Vec2d& p) { return std::abs(p.norm() - r2) < 1e-5; });
 </syntaxhighlight>
-In the code above, we constructed two BallShapes and calculated their difference, giving us the shape of an annulus. We use the annulus shape to obtain the outer domain discretization, and the BallShape with the smaller radius to construct the inner domain discretization. Some care needs to be given to filling the domains, we need to be careful because we want the domains to share the nodes on the boundary $\partial \Omega_2$. Notice that as we discretized the boundary of the annulus we already filled both $\partial \Omega_1$ and $\partial \Omega_2$. Hence, we filter the boundary nodes of the outer domain into the outer boundary and the common boundary, this will prove useful later on.
+In the code above, we constructed two BallShapes and calculated their difference, giving us the shape of an annulus. We use the annulus shape to obtain the outer domain discretization, and the BallShape with the smaller radius to construct the inner domain discretization. Some care needs to be given to filling the domains, because we want the domains to have the same node coordinates on the common boundary $\partial \Omega_2$. Note that as we discretized the boundary of the annulus we already filled both $\partial \Omega_1$ and $\partial \Omega_2$. Hence, we split the boundary nodes of the outer domain into the outer boundary and the common boundary, which will prove useful later on.
 <syntaxhighlight lang="c++" line>
      // Maps node indices of outer domain to their corresponding indices in the inner domain.
@@ Line 44: / Line 44: @@
      inner.fill(fill_engine, h); // this is the inner circle
 </syntaxhighlight>
-The snippet above shows how to add the the nodes on the common boundary to the inner domain. The normals of the shared boundary nodes belonging to the inner domain must point into the opposite direction as the normals belonging to the shared boundary nodes in the outer domain. Additionally a map pointing from outer to inner indexes was made. After this is done insides of both domains are filled. The size of the matrix $M$ of the system must be equal to the sum of points in the inner and outer domains.
+The snippet above shows how to add the the nodes on the common boundary to the inner domain. The normals of the shared boundary nodes belonging to the inner domain must point into the opposite direction as the normals belonging to the shared boundary nodes in the outer domain. Additionally, mapping of node indices from outer domain to the index of the same node in the inner domain was made. After this is done, the insides of both domains are filled. The size of the matrix $M$ of the system must be equal to the sum of numbers of nodes in the inner and outer domains.
 <syntaxhighlight lang="c++" line>
      Eigen::SparseMatrix<double, Eigen::RowMajor> M(N_inner+N_outer, N_inner+N_outer);
      Eigen::VectorXd rhs(N_inner+N_outer);
 </syntaxhighlight>
-Additionaly we must construct two sets of storages, and two sets of operators, one for each domain. The operators of the outer domain must be offset by the number of nodes in the inner domain, this assures that the operators are writing values in the right part of the matrix $M$.
+Additionaly we must construct two sets of storages, and two sets of operators, one for each domain. The operators of the outer domain must be offset by the number of nodes in the inner domain. This assures that the operators are writing values in the correct part of the matrix $M$.
 <syntaxhighlight lang="c++" line>
      auto storage_inner = inner.computeShapes<sh::lap|sh::d1>(wls);
      auto storage_outer = outer.computeShapes<sh::lap|sh::d1>(wls);
@@ Line 70: / Line 69: @@
      M.reserve(reserve);
 </syntaxhighlight>
-Finally we can implement the case, being careful to write in the appropriate parts of the matrix when dealing with shared points. When we filtered the indexes into the common_bnd vector we used the indexes of the nodes in the outer domain, we use the mapping outer_to_inner to obtain the indexes corresponding to the same nodes in the inner domain.
+Finally we can implement the case, being careful to write in the appropriate parts of the matrix when dealing with shared points. When we filtered the indexes into the <code>common_bnd</code> vector we used the indexes of the nodes in the outer domain, we use the mapping <code>outer_to_inner</code> to obtain the indexes corresponding to the same nodes in the inner domain.
 <syntaxhighlight lang="c++" line>
      for (int i : inner.interior()) {
@@ Line 84: / Line 83: @@
          int j = outer_to_inner[i];
          op_inner.value(j) + -1.0*op_outer.value(i, -N_inner+j) = 0.0; // continuity over shared boundary
-         lam1*op_inner.neumann(j, inner.normal(j), N_inner+i) + lam2*op_outer.neumann(i, outer.normal(i)) = 0.0; // flux
+         lam1*op_inner.neumann(j, inner.normal(j), N_inner+i) + lam2*op_outer.neumann(i, outer.normal(i)) = 0.0; // flux continuity
+     }
 </syntaxhighlight>
-The annotated sparse matrix representing the coupled system is shown below
+The annotated sparse matrix representing the coupled system is shown below.
 [[File:couplingMatrix.png|center|600px]]
 The only thing that remains to be done is to solve the system of equations we obtained. Depending on the relative values of $\lambda_1$ and $\lambda_2$ we may obtain two types of solutions

@@ Line 15: / Line 15: @@
 \end{align*}
 We will see that the task is quite simple using Medusa.
 == Discretisation ==
 First we create the domains $\Omega_1$ and $\Omega_2$

@@ Line 1: / Line 1: @@
-In this example we will show how one can solve coupled problems with Medusa library. We will solve the stationary heat equation on concentric circles where the thermal conductivity $\lambda$ is different for each circle. The inner circle will be heated up, while the outside temperature will be held at zero.
+In this example we will show how one can solve coupled problems with the Medusa library. We will solve the stationary heat equation on concentric circles where the thermal conductivity $\lambda$ is different for each circle. The inner circle will be heated up, while the outside temperature will be held at zero.
 [[File:coupledHeat.png|center|300px]]
-In mathematical terms, we want to solve the system
+In mathematical terms, we want to solve the following system for two unknown function $u$ and $v$
 \begin{align*}
      -\lambda_2\Delta v = q_0 \qquad &\text{on} \quad \Omega_2 \\
       \lambda_1\Delta u = 0 \qquad &\text{on} \quad \Omega_1
 \end{align*}
+with boundary conditions
 \begin{align*}
      u = 0 \qquad &\text{on} \quad \partial \Omega_1 \\
      u = v \qquad &\text{on} \quad \partial \Omega_2 \\
-     \lambda_2 \frac{\partial u}{\partial n_1} + \lambda_1 \frac{\partial v}{\partial n_2} = 0 \qquad &\text{on} \quad \partial \Omega_2 \\
+     \lambda_2 \frac{\partial u}{\partial n_1} + \lambda_1 \frac{\partial v}{\partial n_2} = 0 \qquad &\text{on} \quad \partial \Omega_2. \\
 \end{align*}
-We will see that the task is quite simple using Medusa, first we create the domains $\Omega_1$ and $\Omega_2$
+We will see that the task is quite simple using Medusa. First we create the domains $\Omega_1$ and $\Omega_2$
 <syntaxhighlight lang="c++" line>

← Older revision		Revision as of 17:02, 18 April 2019
Line 93:		Line 93:
	</gallery>		</gallery>

−	The example can be found here [~~http~~://~~www~~.~~example~~.~~com link title~~], and the matlab plot script as [~~http~~://~~www~~.~~example~~.~~com link title~~].	+	The example can be found here [https://gitlab.com/e62Lab/medusa/blob/dev/examples/coupled_domains/poisson_coupled_domains.cpp poisson_coupled_domains.cpp], and the matlab plot script as [https://gitlab.com/e62Lab/medusa/blob/dev/examples/coupled_domains/poisson_coupled_domains.m poisson_coupled_domains.m].