Medusa: Coordinate Free Mehless Method implementation - User contributions [en]

File:Maximum value.png

2019-09-17T10:21:37Z

Agrebenc: Agrebenc uploaded a new version of File:Maximum value.png

File uploaded with MsUpload

Heat equation

2019-09-17T09:59:29Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

__FORCETOC__

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).

==Fuction==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t)
\end{align*}
</math>

For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.

<math>
\begin{align*}
relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}}
\end{align*}
</math>

==Weighted Least Squares (WLS)==

In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.

<syntaxhighlight lang="c++" line>
WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);
</syntaxhighlight>

==Finding steady state==

To finish simulation we save values of nodes and compare them after some time steps. Next graph shows us how relative error changes, when we increase number of steps, after which we compare values.

[[File:ending_conditions.png|800px]]

In this graph we can see, that using 8000 time steps is sufficient for the size of our domain, which we will use in further simulations.

==Effect of sigma and weight==

To see the effect of sigma, we ran the same simulation twice and got following graphs:

[[File:different_sigma_1.png|800px]]
[[File:different_sigma_3.png|800px]]

Graphs aren't identical. This is probably due to a slightly different domain. This happens because domain is generated with a bit of randomness.

We can see is that value of sigma has a big effect on relative error.

Weight has simmilar effect on simulation as well.

[[File:different_weights_1.png|800px]]
[[File:different_weights_2.png|800px]]

==Heat transfer==

Now we use sigma 75 and weight 1.0 to solve heat transfer problem. In this case, we don't have analytical sollution. For our boundary condition, we will set up to 0.0, down to 1.0, right to 0.35 and left to 0.75. For the center hole we will use Neumann boundary conditions. We are solving following equations:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial t} = \Delta u(x,y,t)
\end{align*}
</math>

For the domain interior and for domain center we use:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial n} = 0
\end{align*}
</math>

We get the following picture.

[[File:heat_trasfer.png|800px]]

Same simulation with smaller dx value gives us.
[[File:heat_trasfer_2.png|800px]]

We can expect this to have the same error as the simulation with previous function. We can evaluate this by checking what the maximum value is within the domain. Doing so at different numbers of nodes gives us next graph.

[[File:maximum_value.png|800px]]

Maximum value does't exceed 1.0, which it shouldn't.

Go back to [[Medusa#Examples|Examples]].

File:Heat trasfer 2.png

2019-09-17T09:58:31Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

File:Heat trasfer.png

2019-09-17T09:58:13Z

Agrebenc: Agrebenc uploaded a new version of File:Heat trasfer.png

File uploaded with MsUpload

Heat equation

2019-09-17T09:55:27Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

__FORCETOC__

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).

==Fuction==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t)
\end{align*}
</math>

For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.

<math>
\begin{align*}
relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}}
\end{align*}
</math>

==Weighted Least Squares (WLS)==

In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.

<syntaxhighlight lang="c++" line>
WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);
</syntaxhighlight>

==Finding steady state==

To finish simulation we save values of nodes and compare them after some time steps. Next graph shows us how relative error changes, when we increase number of steps, after which we compare values.

[[File:ending_conditions.png|800px]]

In this graph we can see, that using 8000 time steps is sufficient for the size of our domain, which we will use in further simulations.

==Effect of sigma and weight==

To see the effect of sigma, we ran the same simulation twice and got following graphs:

[[File:different_sigma_1.png|800px]]
[[File:different_sigma_3.png|800px]]

Graphs aren't identical. This is probably due to a slightly different domain. This happens because domain is generated with a bit of randomness.

We can see is that value of sigma has a big effect on relative error.

Weight has simmilar effect on simulation as well.

[[File:different_weights_1.png|800px]]
[[File:different_weights_2.png|800px]]

==Heat transfer==

Now we use sigma 75 and weight 1.0 to solve heat transfer problem. In this case, we don't have analytical sollution. For our boundary condition, we will set up to 0.0, down to 1.0, right to 0.35 and left to 0.75. For the center hole we will use Neumann boundary conditions. We are solving following equations:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial t} = \Delta u(x,y,t)
\end{align*}
</math>

For the domain interior and for domain center we use:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial n} = 0
\end{align*}
</math>

We get the following picture.

[[File:heat_trasfer.png|800px]]

We can expect this to have the same error as the simulation with previous function. We can evaluate this by checking what the maximum value is within the domain. Doing so at different numbers of nodes gives us next graph.

[[File:maximum_value.png|800px]]

Maximum value does't exceed 1.0, which it shouldn't.

Go back to [[Medusa#Examples|Examples]].

File:Maximum value.png

2019-09-17T09:55:10Z

Agrebenc: Agrebenc uploaded a new version of File:Maximum value.png

File uploaded with MsUpload

Heat equation

2019-09-17T09:22:49Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

__FORCETOC__

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).

==Fuction==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t)
\end{align*}
</math>

For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.

<math>
\begin{align*}
relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}}
\end{align*}
</math>

==Weighted Least Squares (WLS)==

In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.

<syntaxhighlight lang="c++" line>
WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);
</syntaxhighlight>

==Finding steady state==

To finish simulation we save values of nodes and compare them after some time steps. Next graph shows us how relative error changes, when we increase number of steps, after which we compare values.

[[File:ending_conditions.png|800px]]

In this graph we can see, that using 8000 time steps is sufficient for the size of our domain, which we will use in further simulations.

==Effect of sigma and weight==

To see the effect of sigma, we ran the same simulation twice and got following graphs:

[[File:different_sigma_1.png|800px]]
[[File:different_sigma_3.png|800px]]

Graphs aren't identical. This is probably due to a slightly different domain. This happens because domain is generated with a bit of randomness.

We can see is that value of sigma has a big effect on relative error.

Weight has simmilar effect on simulation as well.

[[File:different_weights_1.png|800px]]
[[File:different_weights_2.png|800px]]

==Heat transfer==

Now we use sigma 75 and weight 1.0 to solve heat transfer problem. In this case, we don't have analytical sollution. For our boundary condition, we will set up to 0.0, down to 1.0, right to 0.35 and left to 0.75. For the center hole we will use Neumann boundary conditions. We are solving following equations:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial t} = \Delta u(x,y,t)
\end{align*}
</math>

For the domain interior and for domain center we use:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial n} = 0
\end{align*}
</math>

We get the following picture.

[[File:heat_trasfer.png|800px]]

We can expect this to have the same error as the simulation with previous function. We can evaluate this by checking what the maximum value is within the domain. Doing so at different numbers of nodes gives us next graph.

[[File:maximum_value.png|800px]]

Maximum value does't exceed 1.0.

Go back to [[Medusa#Examples|Examples]].

File:Different weights 2.png

2019-09-17T09:05:30Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

File:Different weights 1.png

2019-09-17T09:05:22Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

Heat equation

2019-09-17T08:42:47Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

__FORCETOC__

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).

==Fuction==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t)
\end{align*}
</math>

For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.

<math>
\begin{align*}
relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}}
\end{align*}
</math>

==Weighted Least Squares (WLS)==

In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.

<syntaxhighlight lang="c++" line>
WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);
</syntaxhighlight>

==Finding steady state==

To finish simulation we save values of nodes and compare them after some time steps. Next graph shows us how relative error changes, when we increase number of steps, after which we compare values.

[[File:ending_conditions.png|800px]]

In this graph we can see, that using 8000 time steps is sufficient for the size of our domain, which we will use in further simulations.

==Effect of sigma and weight==

To see the effect of sigma, we ran the same simulation twice and got following graphs:

[[File:different_sigma_1.png|800px]]
[[File:different_sigma_3.png|800px]]

(Graphs aren't identical. This is probably due to a slightly different domain. This happens becouse domain is generated with a bit of randomness.)?

We can see is that value of sigma has a big effect on relative error.

Weight - TODO.

==Heat transfer==

Now we use sigma 75 and weight 1.0 to solve heat transfer problem. In this case, we don't have analytical sollution. For our boundary condition, we will set up to 0.0, down to 1.0, right to 0.35 and left to 0.75. For the center hole we will use Neumann boundary conditions. We are solving following equations:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial t} = \Delta u(x,y,t)
\end{align*}
</math>

For the domain interior and for domain center we use:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial n} = 0
\end{align*}
</math>

We get the following picture.

[[File:heat_trasfer.png|800px]]

We can expect this to have the same error as the simulation with previous function. We can evaluate this by checking what the maximum value is within the domain. Doing so at different numbers of nodes gives us next graph.

[[File:maximum_value.png|800px]]

Maximum value does't exceed 1.0.

Go back to [[Medusa#Examples|Examples]].

File:Different sigma 3.png

2019-09-17T08:42:41Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

File:Different sigma 2.png

2019-09-17T08:41:11Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 2.png

File uploaded with MsUpload

Heat equation

2019-09-17T08:40:04Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

__FORCETOC__

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).

==Fuction==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t)
\end{align*}
</math>

For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.

<math>
\begin{align*}
relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}}
\end{align*}
</math>

==Weighted Least Squares (WLS)==

In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.

<syntaxhighlight lang="c++" line>
WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);
</syntaxhighlight>

==Finding steady state==

To finish simulation we save values of nodes and compare them after some time steps. Next graph shows us how relative error changes, when we increase number of steps, after which we compare values.

[[File:ending_conditions.png|800px]]

In this graph we can see, that using 8000 time steps is sufficient for the size of our domain, which we will use in further simulations.

==Effect of sigma and weight==

To see the effect of sigma, we ran the same simulation twice and got following graphs:

[[File:different_sigma_1.png|800px]]
[[File:different_sigma_2.png|800px]]

(Graphs aren't identical. This is probably due to a slightly different domain. This happens becouse domain is generated with a bit of randomness.)?

We can see is that value of sigma has a big effect on relative error.

Weight - TODO.

==Heat transfer==

Now we use sigma 75 and weight 1.0 to solve heat transfer problem. In this case, we don't have analytical sollution. For our boundary condition, we will set up to 0.0, down to 1.0, right to 0.35 and left to 0.75. For the center hole we will use Neumann boundary conditions. We are solving following equations:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial t} = \Delta u(x,y,t)
\end{align*}
</math>

For the domain interior and for domain center we use:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial n} = 0
\end{align*}
</math>

We get the following picture.

[[File:heat_trasfer.png|800px]]

We can expect this to have the same error as the simulation with previous function. We can evaluate this by checking what the maximum value is within the domain. Doing so at different numbers of nodes gives us next graph.

[[File:maximum_value.png|800px]]

Maximum value does't exceed 1.0.

Go back to [[Medusa#Examples|Examples]].

File:Different sigma 1.png

2019-09-17T08:39:12Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 1.png

File uploaded with MsUpload

File:Different sigma 2.png

2019-09-17T08:39:03Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 2.png

File uploaded with MsUpload

File:Different sigma 2.png

2019-09-17T08:37:14Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 2.png

File uploaded with MsUpload

File:Different sigma 1.png

2019-09-17T08:37:03Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 1.png

File uploaded with MsUpload

File:Different sigma 2.png

2019-09-17T08:35:47Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 2.png

File uploaded with MsUpload

File:Different sigma 1.png

2019-09-17T08:35:38Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 1.png

File uploaded with MsUpload

File:Different sigma 2.png

2019-09-17T07:10:02Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 2.png

File uploaded with MsUpload

File:Different sigma 1.png

2019-09-17T07:09:54Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 1.png

File uploaded with MsUpload

File:Different sigma 1.png

2019-09-17T07:09:43Z

Agrebenc: Agrebenc uploaded a new version of File:Different sigma 1.png

File uploaded with MsUpload

Heat equation

2019-09-13T11:29:05Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

__FORCETOC__

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).

==Fuction==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t)
\end{align*}
</math>

For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.

<math>
\begin{align*}
relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}}
\end{align*}
</math>

==Weighted Least Squares (WLS)==

In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.

<syntaxhighlight lang="c++" line>
WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);
</syntaxhighlight>

==Finding steady state==

To finish simulation we save values of nodes and compare them after some time steps. Next graph shows us how relative error changes, when we increase number of steps, after which we compare values.

[[File:ending_conditions.png|800px]]

In this graph we can see, that using 8000 time steps is sufficient for the size of our domain, which we will use in further simulations.

==Effect of sigma and weight==

To see the effect of sigma, we ran the same simulation twice and got following graphs:

[[File:different_sigma_1.png|800px]] [[File:different_sigma_2.png|800px]]

(Graphs aren't identical. This is probably due to a slightly different domain. This happens becouse domain is generated with a bit of randomness.)?

We can see is that value of sigma has a big effect on relative error.

Weight - TODO.

==Heat transfer==

Now we use sigma 75 and weight 1.0 to solve heat transfer problem. In this case, we don't have analytical sollution. For our boundary condition, we will set up to 0.0, down to 1.0, right to 0.35 and left to 0.75. For the center hole we will use Neumann boundary conditions. We are solving following equations:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial t} = \Delta u(x,y,t)
\end{align*}
</math>

For the domain interior and for domain center we use:

<math>
\begin{align*}
\frac{\partial u (x,y,t)}{\partial n} = 0
\end{align*}
</math>

We get the following picture.

[[File:heat_trasfer.png|800px]]

We can expect this to have the same error as the simulation with previous function. We can evaluate this by checking what the maximum value is within the domain. Doing so at different numbers of nodes gives us next graph.

[[File:maximum_value.png|800px]]

Maximum value does't exceed 1.0.

Go back to [[Medusa#Examples|Examples]].

File:Maximum value.png

2019-09-13T09:50:44Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

File:Heat trasfer.png

2019-09-13T09:45:36Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

Heat equation

2019-09-12T09:23:57Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

__FORCETOC__

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).

==Fuction==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t)
\end{align*}
</math>

For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.

<math>
\begin{align*}
relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}}
\end{align*}
</math>

==Weighted Least Squares (WLS)==

In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.

<syntaxhighlight lang="c++" line>
WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);
</syntaxhighlight>

==Finding steady state==

To finish simulation we save values of nodes and compare them after some time steps. Next graph shows us how relative error changes, when we increase number of steps, after which we compare values.

[[File:ending_conditions.png|800px]]

In this graph we can see, that using 8000 time steps is sufficient for the size of our domain, which we will use in further simulations.

==Effect of sigma and weight==

To see the effect of sigma, we ran the same simulation twice and got following graphs:

[[File:different_sigma_1.png|800px]] [[File:different_sigma_2.png|800px]]

(Graphs aren't identical. This is probably due to a slightly different domain. This happens becouse domain is generated with a bit of randomness.)?

We can see is that value of sigma has a big effect on relative error.

Go back to [[Medusa#Examples|Examples]].

Heat equation

2019-09-12T09:23:37Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

This example shows how explicit method of solving partial differential equations is used within the Medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

__FORCETOC__

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them. Domain boundary is separated in to 5 areas: up, down, left, right and center. The idea is to have Dirichlet boundary condition on the square (up, down, left and right) and Neumann boundary conditions in the middle (center).

==Fuction==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to $u\_analytic$. With this we also calculate Dirichlet boundary conditions. Then we set domain interior and nodes with Neumann bondary conditions to arbitrary value, in this case 0. We proceed to explicitly solve the domain. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t)
\end{align*}
</math>

For Neumann bondary conditions we set $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

After the simulation is finished, we can easily calculate relative error. Values of nodes are saved in a vector $u1$, so we can use norm of $|u1 - u\_analytic|$ to calculate error, then divide it by norm of $u1$ to calculate relative error.

<math>
\begin{align*}
relative\_error = \frac{\sqrt{\sum{|u1-u\_analytic|^2}}}{\sqrt{\sum{|u1|^2}}}
\end{align*}
</math>

==Weighted Least Squares (WLS)==

In this simulation, we use Gaussians for our WLS. Because of scattred nodes in our domain, we have to use 9 basis. The support size we chose is 27. For weight we have chosen 1.0 and for Gaussian weight (sigma) 75. Later we will look how change of sigma and weight effects relative error.

<syntaxhighlight lang="c++" line>
WLS<Gaussians<Vec2d>, GaussianWeight<Vec2d>, ScaleToClosest> wls({9, 75.0}, 1.0);
</syntaxhighlight>

==Finding steady state==

To finish simulation we save values of nodes and compare them after some time steps. Next graph shows us how relative error changes, when we increase number of steps, after which we compare values.

[[File:ending_conditions.png|800px]]

In this graph we can see, that using 8000 time steps is sufficient for the size of our domain, which we will use in further simulations.

==Effect of sigma and weight==

To see the effect of sigma, we ran the same simulation twice and got following graphs:

[[File:different_sigma_1.png|800px]] [[File:different_sigma_2.png|800px]]

(Graphs aren't identical. This is probably due to a slightly different domain. This happens becouse domain is generated with a bit of randomness.)?
We can see is that value of sigma has a big effect on relative error.

Go back to [[Medusa#Examples|Examples]].

File:Different sigma 2.png

2019-09-12T09:19:32Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

File:Different sigma 1.png

2019-09-12T09:19:19Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

Heat equation

2019-09-12T08:37:11Z

Agrebenc:

File:Ending conditions.png

2019-09-12T08:06:26Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

Heat equation

2019-08-19T13:16:12Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

__FORCETOC__

==Purpose of Heat equation example==

This example shows how explicit method of solving partial differential equations is used within the medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them.

<syntaxhighlight lang="c++" line>
BoxShape<Vec2d> box(-3.0, 3.0);
DomainDiscretization<Vec2d> domain=box.discretizeBoundaryWithStep(dx);
BallShape<Vec2d> ball(0,1.0);
DomainDiscretization<Vec2d> domain_empty=ball.discretizeBoundaryWithStep(dx,-5);
domain -=domain_empty;
GeneralFill<Vec2d> fill_engine;
fill_engine(domain, dx);
</syntaxhighlight>

Domain boundary is separated in to 5 areas: up, down, left, right and center.

<syntaxhighlight lang="c++" line>
auto internal = domain.interior();
auto up = domain.types().filter([](int i) { return i == -4 ; });
auto down = domain.types().filter([](int i) { return i == -3 ; });
auto right = domain.types().filter([](int i) { return i == -2; });
auto left = domain.types().filter([](int i) { return i == -1; });
auto center=domain.types().filter([](int i){ return i == -5;});
</syntaxhighlight>

The idea of this domain is to have Dirichlet boundary condition on the square. Furthermore, they are constant over the course of simulation. In the center, the boundary conditions are Neumann.

==Fuction, solution and comparison==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to u1_analytic. With this we also calculate Dirichlet boundary conditions. Then we set interior and center of domain to arbitrary value, in this case 0. Then we proceed to explicitly solve the domain interior and center. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*\bigg(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y)\bigg)+u(x,y,t)
\end{align*}
</math>

For domain center we set value of $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

While doing this, we also compare previous values u1 to new ones u2. If there are differences bigger than dt*1E-5, we continue the simulation. Otherwise, we stop it and compare values in nodes to analytical ones.

<syntaxhighlight lang="c++" line>
while (ending){
if (iteracije%izpis==0) {
std::cout << iteracije << std::endl;
}
iteracije++;
#pragma omp parallel for
for (int j = 0; j < internal.size(); ++j) {
int i = internal[j];
u2[i] = dt *( op.lap(u1, i) +2*PI*PI*function(domain.pos(i,0),domain.pos(i, 1)))+ u1[i];
}
#pragma omp parallel for
for (int j = 0; j < center.size(); ++j) {
int i = center[j];
auto vec=domain.normal(i);
double x=domain.pos(i,0);
double y=domain.pos(i, 1);
u2[i]=op.neumann(u1,i,vec,vec[0]*PI*cos(PI*x)*sin(PI*y)+vec[1]*PI*cos(PI*y)*sin(PI*x));
}
for (int i:internal+center){
double razlika=abs(u1[i]-u2[i]);
if (razlika>konec){
break;
}
ending = false;
}
u1.swap(u2);
}
</syntaxhighlight>

We ran simulation a couple of times using different dx values. Changing dx value changed number of total nodes as well. We made logarithmic graph of relative error per number of nodes.

[[File:error_total_nodes.png|800px]]

Go back to [[Medusa#Examples|Examples]].

File:Error total nodes.png

2019-08-19T12:43:44Z

Agrebenc: File uploaded with MsUpload

File uploaded with MsUpload

Heat equation

2019-08-19T12:42:39Z

Agrebenc:

Heat equation

2019-08-19T12:21:54Z

Agrebenc:

Heat equation

2019-08-19T11:07:57Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

__FORCETOC__

==Purpose of Heat equation example==

This example shows how explicit method of solving partial differential equations is used within the medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them.

<syntaxhighlight lang="c++" line>
BoxShape<Vec2d> box(-3.0, 3.0);
DomainDiscretization<Vec2d> domain=box.discretizeBoundaryWithStep(dx);
BallShape<Vec2d> ball(0,1.0);
DomainDiscretization<Vec2d> domain_empty=ball.discretizeBoundaryWithStep(dx,-5);
domain -=domain_empty;
GeneralFill<Vec2d> fill_engine;
fill_engine(domain, dx);
</syntaxhighlight>

Domain boundary is separated in to 5 areas: up, down, left, right and center.

<syntaxhighlight lang="c++" line>
auto internal = domain.interior();
auto up = domain.types().filter([](int i) { return i == -4 ; });
auto down = domain.types().filter([](int i) { return i == -3 ; });
auto right = domain.types().filter([](int i) { return i == -2; });
auto left = domain.types().filter([](int i) { return i == -1; });
auto center=domain.types().filter([](int i){ return i == -5;});
</syntaxhighlight>

The idea of this domain is to have Dirichlet boundary condition on the square. Furthermore, they are constant over the course of simulation. In the center, the boundary conditions are Neumann.

==Fuction, solution and comparison==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to u1_analytic. With this we also calculate Dirichlet boundary conditions. Then we set interior and center of domain to arbitrary value, in this case 0. Then we proceed to explicitly solve the domain interior and center. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y))+u(x,y,t)
\end{align*}
</math>

For domain center we set value of $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

While doing this, we also compare previous values u1 to new ones u2. If there are differences bigger than dt*1E-5, we continue the simulation. Otherwise, we stop it and compare values in nodes to analytical ones.

<syntaxhighlight lang="c++" line>
while (ending){
if (iteracije%izpis==0) {
std::cout << iteracije << std::endl;
}
iteracije++;
#pragma omp parallel for
for (int j = 0; j < internal.size(); ++j) {
int i = internal[j];
u2[i] = dt *( op.lap(u1, i) +2*PI*PI*function(domain.pos(i,0),domain.pos(i, 1)))+ u1[i];
}
#pragma omp parallel for
for (int j = 0; j < center.size(); ++j) {
int i = center[j];
auto vec=domain.normal(i);
double x=domain.pos(i,0);
double y=domain.pos(i, 1);
u2[i]=op.neumann(u1,i,vec,vec[0]*PI*cos(PI*x)*sin(PI*y)+vec[1]*PI*cos(PI*y)*sin(PI*x));
}
for (int i:internal+center){
double razlika=abs(u1[i]-u2[i]);
if (razlika>konec){
break;
}
ending = false;
}
u1.swap(u2);
}
</syntaxhighlight>

Go back to [[Medusa#Examples|Examples]].

Heat equation

2019-08-19T11:07:40Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].
Page in progress!
__FORCETOC__

==Purpose of Heat equation example==

This example shows how explicit method of solving partial differential equations is used within the medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them.

<syntaxhighlight lang="c++" line>
BoxShape<Vec2d> box(-3.0, 3.0);
DomainDiscretization<Vec2d> domain=box.discretizeBoundaryWithStep(dx);
BallShape<Vec2d> ball(0,1.0);
DomainDiscretization<Vec2d> domain_empty=ball.discretizeBoundaryWithStep(dx,-5);
domain -=domain_empty;
GeneralFill<Vec2d> fill_engine;
fill_engine(domain, dx);
</syntaxhighlight>

Domain boundary is separated in to 5 areas: up, down, left, right and center.

<syntaxhighlight lang="c++" line>
auto internal = domain.interior();
auto up = domain.types().filter([](int i) { return i == -4 ; });
auto down = domain.types().filter([](int i) { return i == -3 ; });
auto right = domain.types().filter([](int i) { return i == -2; });
auto left = domain.types().filter([](int i) { return i == -1; });
auto center=domain.types().filter([](int i){ return i == -5;});
</syntaxhighlight>

The idea of this domain is to have Dirichlet boundary condition on the square. Furthermore, they are constant over the course of simulation. In the center, the boundary conditions are Neumann.

==Fuction, solution and comparison==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to u1_analytic. With this we also calculate Dirichlet boundary conditions. Then we set interior and center of domain to arbitrary value, in this case 0. Then we proceed to explicitly solve the domain interior and center. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y))+u(x,y,t)
\end{align*}
</math>

For domain center we set value of $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

While doing this, we also compare previous values u1 to new ones u2. If there are differences bigger than dt*1E-5, we continue the simulation. Otherwise, we stop it and compare values in nodes to analytical ones.

<syntaxhighlight lang="c++" line>
while (ending){
if (iteracije%izpis==0) {
std::cout << iteracije << std::endl;
}
iteracije++;
#pragma omp parallel for
for (int j = 0; j < internal.size(); ++j) {
int i = internal[j];
u2[i] = dt *( op.lap(u1, i) +2*PI*PI*function(domain.pos(i,0),domain.pos(i, 1)))+ u1[i];
}
#pragma omp parallel for
for (int j = 0; j < center.size(); ++j) {
int i = center[j];
auto vec=domain.normal(i);
double x=domain.pos(i,0);
double y=domain.pos(i, 1);
u2[i]=op.neumann(u1,i,vec,vec[0]*PI*cos(PI*x)*sin(PI*y)+vec[1]*PI*cos(PI*y)*sin(PI*x));
}
for (int i:internal+center){
double razlika=abs(u1[i]-u2[i]);
if (razlika>konec){
break;
}
ending = false;
}
u1.swap(u2);
}
</syntaxhighlight>

Go back to [[Medusa#Examples|Examples]].

Heat equation

2019-08-19T10:34:22Z

Agrebenc:

Go back to [[Medusa#Examples|Examples]].

__FORCETOC__

==Purpose of Heat equation example==

This example shows how explicit method of solving partial differential equations is used within the medusa library. Firstly, this example compares explicit solutions to analytical, then it shows how it can be used to solve heat transfer equation.

==Domain==

In this example, a square domain with a cut-out hole in the middle is used. The hole has radius equal to 1/3 of length of a square. Domain nodes are placed randomly with about dx length distance between them.

<syntaxhighlight lang="c++" line>
BoxShape<Vec2d> box(-3.0, 3.0);
DomainDiscretization<Vec2d> domain=box.discretizeBoundaryWithStep(dx);
BallShape<Vec2d> ball(0,1.0);
DomainDiscretization<Vec2d> domain_empty=ball.discretizeBoundaryWithStep(dx,-5);
domain -=domain_empty;
GeneralFill<Vec2d> fill_engine;
fill_engine(domain, dx);
</syntaxhighlight>

Domain boundary is separated in to 5 areas: up, down, left, right and center.

<syntaxhighlight lang="c++" line>
auto internal = domain.interior();
auto up = domain.types().filter([](int i) { return i == -4 ; });
auto down = domain.types().filter([](int i) { return i == -3 ; });
auto right = domain.types().filter([](int i) { return i == -2; });
auto left = domain.types().filter([](int i) { return i == -1; });
auto center=domain.types().filter([](int i){ return i == -5;});
</syntaxhighlight>

The idea of this domain is to have Dirichlet boundary condition on the square. Furthermore, they are constant over the course of simulation. In the center, the boundary conditions are Neumann.

==Fuction for comparison of solutions==

The function used for this example is $u(x,y) = \sin(\pi x)\sin(\pi y)$. Laplace operator of this function is $\Delta u(x,y)=-2\pi^2\sin(\pi x)\sin(\pi y)$.
For comparing explicit solution to analytical, we solve function for each node and save solution to u1_analytic. With this we also calculate Dirichlet boundary conditions. Then we set interior and center of domain to arbitrary value, in this case 0. Then we proceed to explicitly solve the domain interior and center. In the domain interior we are solving the following equation:

<math>
\begin{align*}
u (x,y,t+1) = dt*(\Delta u(x,y,t) + 2\pi^2\sin(\pi x)\sin(\pi y))+u(x,y,t)
\end{align*}
</math>

For domain center we set value of $u(x,y,t+1)$ in a way that it satisfies the following equation:

<math>
\frac{\partial u(x,y,t+1)}{\partial n}=n_x \pi \cos(\pi x) \sin(\pi y)+n_y \pi \sin(\pi x) \cos(\pi y)
</math>

While doing this, we also compare previous values u1 to new ones u2. If there are differences bigger than dt*1E-5, we continue the simulation. Otherwise, we stop it and compare values in nodes to analytical ones.

<syntaxhighlight lang="c++" line>
while (ending){
if (iteracije%izpis==0) {
std::cout << iteracije << std::endl;
}
iteracije++;
#pragma omp parallel for
for (int j = 0; j < internal.size(); ++j) {
int i = internal[j];
u2[i] = dt *( op.lap(u1, i) +2*PI*PI*function(domain.pos(i,0),domain.pos(i, 1)))+ u1[i];
}
#pragma omp parallel for
for (int j = 0; j < center.size(); ++j) {
int i = center[j];
auto vec=domain.normal(i);
double x=domain.pos(i,0);
double y=domain.pos(i, 1);
u2[i]=op.neumann(u1,i,vec,vec[0]*PI*cos(PI*x)*sin(PI*y)+vec[1]*PI*cos(PI*y)*sin(PI*x));
}
for (int i:internal+center){
double razlika=abs(u1[i]-u2[i]);
if (razlika>konec){
break;
}
ending = false;
}
u1.swap(u2);
}
</syntaxhighlight>

Go back to [[Medusa#Examples|Examples]].

Heat equation

2019-08-19T09:55:05Z

Agrebenc: Created page with "Go back to Examples. __FORCETOC__ ==Purpose of Heat equation example== This example shows how explicit method of solving partial differential equations..."