StaticVectorSolver::Solver1< dim, stage > Class Template Reference

Solves static vector boundary value problem. More...

#include <static_vector_solver_i.hpp>

Inheritance diagram for StaticVectorSolver::Solver1< dim, stage >:

Public Member Functions
	Solver1 (unsigned int p, unsigned int mapping_degree, unsigned int type_of_pde_rhs=3, double eta_squared=0.0, std::string fname="data", const Function< dim > *exact_solution=nullptr, bool print_time_tables=false, bool project_exact_solution=false, bool write_higher_order_cells=false)
	The only constructor. More...
virtual void	make_mesh ()=0
	Initializes the data member StaticVectorSolver::Solver1::triangulation. More...
virtual void	fill_dirichlet_stack ()=0
	Initializes the data member StaticVectorSolver::Solver1::dirichlet_stack. More...
virtual void	solve ()=0
	Solves the system of linear equations.
void	setup ()
	Initializes system matrix and the right-hand side vector. More...
void	assemble ()
	Assembles the system matrix and the right-hand side vector.
void	compute_error_norms ()
	Computes error norms.
void	project_exact_solution_fcn ()
	Projects exact solution. More...
void	save () const
	Saves simulation results into a vtk or vtu file. More...
void	save_matrix_and_rhs_to_csv (std::string fname) const
	Saves the system matrix and the right-hand side into a csv file. More...
void	clear ()
	Releases computer memory associated with system matrix and right-hand side.
const Triangulation< dim > &	get_tria () const
	Returns a reference to triangulation.
const DoFHandler< dim > &	get_dof_handler () const
	Returns a reference to dof handler.
const Vector< double > &	get_solution () const
	Returns a reference to solution.
unsigned int	get_n_cells () const
	Returns the number of active cells in the mesh.
unsigned int	get_n_dofs () const
	Returns the total amount of the degrees of freedom.
unsigned int	get_rhs_type () const
	Returns the value of type_of_pde_rhs.
double	get_L2_norm () const
	Returns \(L^2\) error norm.
double	get_Linfty_norm () const
	Returns \(L^{\infty}\) error norm.
unsigned int	get_mapping_degree () const
	Returns degree of the interpolating Lagrange polynomials used for mapping from the reference cell to the real mesh cell and back.
void	run ()
	Runs the simulation. More...

Protected Attributes
std::map< types::boundary_id, const Function< dim > * >	dirichlet_stack
	A map that contains pairs of boundary IDs and the corresponding Dirichlet boundary conditions. All boundary IDs must be odd.
Triangulation< dim >	triangulation
	The mesh.
const FE_Nedelec< dim >	fe
	The finite elements.
DoFHandler< dim >	dof_handler
	The dof handler.
Vector< double >	solution
	The solution vector, that is, degrees of freedom yielded by the simulation.
Vector< double >	projected_exact_solution
	The projected exact solution vector.
AffineConstraints< double >	constraints
	The constraints associated with the Dirichlet boundary conditions.
SparsityPattern	sparsity_pattern
	The sparsity pattern of the system matrix.
SparseMatrix< double >	system_matrix
	The system matrix.
Vector< double >	system_rhs
	The system right-hand side vector.
double	L2_norm
	The \(L^2\) norm.
double	Linfty_norm
	The \(L^{\infty}\) norm.

Detailed Description

template<int dim, int stage = 1>
class StaticVectorSolver::Solver1< dim, stage >

Solves static vector boundary value problem.

Implements the following recipes:

• (1) Recipe for static vector solver in 3D
• (2) Recipe for static vector solver in 2D
• (3) Recipe for static vector solver in 3D (current vect. potential)

This class template is intended to be a general solver for problems in magnetostatics that can be formulated in terms of the magnetic vector potential, \(\vec{A}\). It can also be used for calculating the current vector potential, \(\vec{T}\), i.e., converting a closed-form analytical expression for \(\vec{J}_f\) into \(\vec{T}\) expressed as a finite-element field function. Such calculated \(\vec{T}\) can be used as an input for StaticVectorSolver::Solver2. The Bossavit's diagrams below illustrate the partial differential equations that can be solved with a help of this class template. In all five cases the vector potential is modeled by the FE_Nedelec finite elements.

A user of this class is supposed to do the following.

• Derive a class from StaticVectorSolver::Solver1.
• Override the virtual destructor.
• Call constructor StaticVectorSolver::Solver1 and pass to it the degree of interpolating polynomials and other arguments.
• Override member functions

  virtual void make_mesh() = 0;
  virtual void fill_dirichlet_stack() = 0;
  virtual void solve() = 0;

• Implement the member functions

  void value_list(...);

  of the classes StaticVectorSolver::TheCoefficient, StaticVectorSolver::PdeRhs, StaticVectorSolver::Gamma, StaticVectorSolver::RobinRhs, and StaticVectorSolver::FreeSurfaceCurrent.
• Implement member function

  double value(...) const;

  of the class StaticVectorSolver::Weight.
• Explicitly call

  StaticVectorSolver::run();

  Alternatively, the user may call individual member functions (such as make_mesh(), fill_dirichlet_stack(), setup(), etc.) in a proper order.

The boundaries of the mesh must be labeled such that the boundary_id() member function of a face object returns the corresponding boundary ID.

The boundary ID's must obey the following convention.

• The Dirichlet boundary conditions are applied on the boundaries with odd boundary ID's.
• The Robin boundary conditions are applied on the boundaries with even boundary ID's. The boundary ID's in this case must be greater than zero.
• No boundary condition is applied on a boundary with zero ID. Applying no boundary condition is as good as applying the homogeneous Neumann boundary condition, \(\big(1/\mu\big)\hat{n}\times\big(\vec{\nabla}\times\vec{A}\big)=0\), as implicitly implied by the first term of the functional. This boundary condition can also be imposed by assigning to a boundary an even ID greater than zero, and setting \( \gamma \) and \(\vec{Q}\) to zero in the value_list methods of the classes StaticVectorSolver::Gamma and StaticVectorSolver::RobinRhs.

When solving for the magnetic vector potential, \(\vec{A}\), the following modes of operation are available:

• type_of_pde_rhs = 0. There is no volume free-current density in the problem domain. The surface free-current density, \(\vec{K}_f\), can be present on interfaces. The right-hand side of the partial differential equation equals zero, i.e.,

  \[ \vec{\nabla}\times\bigg(\dfrac{1}{\mu_0} \vec{\nabla}\times\vec{A}\bigg) + \eta^2 \vec{A} = 0 \]

  in a three-dimensional space and

  \[ \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu_0} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg)+\eta^2\vec{A}=0 \]

  in a two-dimensional space. The data provided by StaticVectorSolver::PdeRhs is not used. This mode allows saving simulation time on computing the integrals associated with \(\vec{J}_f\).
• type_of_pde_rhs = 1. The data provided by StaticVectorSolver::PdeRhs is interpreted as the free-current density, \(\vec{J}_f\), i.e.,

  \[\vec{\nabla} \times\bigg(\dfrac{1}{\mu_0}\vec{\nabla}\times\vec{A}\bigg) +\eta^2 \vec{A} = \vec{J}_f \]

  in a three-dimensional space and

  \[ \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu_0} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg)+ \eta^2 \vec{A} = \vec{J}_f \]

  in a two-dimensional space. The corresponding therm of the functional is

  \[ \iiint_{\Omega} \vec{J}_f \cdot \vec{A} dV \]

  in a three-dimensional space and

  \[ \iint_{\Omega} \vec{J}_f \cdot \vec{A} dS \]

  in a two-dimensional space.
• type_of_pde_rhs = 2. The data provided by StaticVectorSolver::PdeRhs is interpreted as vector current potential, \(\vec{T}\), i.e.,

  \[\vec{\nabla} \times\bigg(\dfrac{1}{\mu_0}\vec{\nabla}\times\vec{A}\bigg) +\eta^2 \vec{A} = \vec{\nabla}\times\vec{T} \]

  in a three-dimensional space and

  \[ \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu_0} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg)+ \eta^2 \vec{A} = \vec{\nabla}\overset{V}{\times} T \]

  in a two-dimensional space. Then the corresponding term of the functional is

  \[ \iiint_{\Omega} \vec{T} \cdot \bigg(\vec{\nabla}\times\vec{A}\bigg) dV \]

  in a three-dimensional space and

  \[ \iint_{\Omega} T \bigg(\vec{\nabla}\overset{S}{\times}\vec{A}\bigg) dS \]

  in a two-dimensional space.
• type_of_pde_rhs = 3. The data provided by StaticVectorSolver::PdeRhs is interpreted as vector current potential, \(\vec{T}\), i.e.,

  \[ \vec{\nabla}\times\bigg(\dfrac{1}{\mu_0}\vec{\nabla}\times\vec{A}\bigg) +\eta^2 \vec{A}=\vec{\nabla}\times\vec{T} \]

  in a three-dimensional space and

  \[ \vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu_0} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg)+\eta^2\vec{A}= \vec{\nabla}\overset{V}{\times} T\]

  in a two-dimensional space. Then the corresponding terms of the functional are

  \[\iiint_{\Omega} \vec{T}\cdot\bigg(\vec{\nabla}\times\vec{A}\bigg) dV - \underbrace{ \iint_{\Gamma_{\Omega}}\vec{T}\cdot\bigg(\hat{n}\times\vec{A}\bigg)dS }_{\text{Boundary integral}} \]

  in a three dimensional space and

  \[ \iint_{\Omega}T\bigg(\vec{\nabla}\overset{S}{\times}\vec{A}\bigg)dS- \underbrace{ \int_{\Gamma_{\Omega}}T\bigg(\hat{n}\overset{S}{\times}\vec{A}\bigg)dl }_{\text{Boundary integral}} \]

  in a two-dimensional space.

The mode type_of_pde_rhs = 2 differs from the mode type_of_pde_rhs = 3 at one point only: the boundary integral is not calculated if type_of_pde_rhs = 2. This can help to reduce simulation time a bit if the current vector potential, \(\vec{T}\), is set to zero by the homogeneous Dirichlet boundary condition.

The same four modes are available when solving for the current vector potential, \(\vec{T}\). In the case of \(\vec{T}\), however, the first two modes, i.e., type_of_pde_rhs=0 and type_of_pde_rhs=1, do not make much sense: the curl of the free-current density must be present on the right-hand side of the partial differential equation. Recall that the current vector potential in two-dimensional problems, \(T\), is always an out-of plane vector, i.e., a scalar. It is described by the div-grad equation, not by the curl-curl equation. That is to say, one must use StaticScalarSolver::Solver to solve for the two-dimensional current vector potential. It has absolutely nothing to do with the StaticVectorSolver::Solver1 described on this page. When solving for the three-dimensional current vector potential, \(\vec{T}\), the following two modes make sense:

• type_of_pde_rhs = 2. The data provided by StaticVectorSolver::PdeRhs is interpreted as free-current density, \(\vec{J}_f\), i.e.,

  \[ \vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) + \eta^2 \vec{T} = \vec{\nabla}\times\vec{J}_f \]

  in a three-dimensional space. Then the corresponding term of the functional is

  \[ \iiint_{\Omega} \vec{J}_f\cdot\bigg(\vec{\nabla}\times\vec{T}\bigg) dV. \]

• type_of_pde_rhs = 3. The data provided by StaticVectorSolver::PdeRhs is interpreted as free-current density, \(\vec{J}_f\),

  \[ \vec{\nabla}\times\bigg(\vec{\nabla}\times\vec{T}\bigg) +\eta^2 \vec{T} = \vec{\nabla}\times\vec{J}_f \]

  in a three-dimensional space. Then the corresponding terms of the functional are

  \[ \iiint_{\Omega} \vec{J}_f\cdot\bigg(\vec{\nabla}\times\vec{T}\bigg) dV - \underbrace{ \iint_{\Gamma_{\Omega}}\vec{J}_f\cdot\bigg(\hat{n}\times\vec{T}\bigg)dS. }_{\text{Boundary integral}} \]

Here again the mode type_of_pde_rhs = 3 differs from the mode type_of_pde_rhs = 2 by the boundary integral in the functional.

Note, the code that implements the solver for \(\vec{A}\) is identical to the code that implements the solver for \(\vec{T}\) as the list of arguments of the constructor of this class template does not contain an argument that toggles between two modes: "solving for A" mode and "solving for T" mode. The user toggles between these two modes by feeding the right data through StaticVectorSolver::PdeRhs and StaticVectorSolver::TheCoefficient and by not using the modes type_of_pde_rhs = 0 and type_of_pde_rhs = 1 when solving for \(\vec{T}\). If the magnetic vector potential, \(\vec{A}\), is being computed, StaticVectorSolver::PdeRhs::value_list must return the values of \(\vec{T}\) (type_of_pde_rhs = 2, type_of_pde_rhs = 3) or values of \(\vec{J}_f\) (type_of_pde_rhs 1). If the current vector potential, \(\vec{T}\), is being computed, dim must equal 3 and StaticVectorSolver::PdeRhs::value_list must return the values of \(\vec{J}_f\) (type_of_pde_rhs = 2, type_of_pde_rhs = 3). The StaticVectorSolver::TheCoefficient must return 1.0 when solving for \(\vec{T}\) and \(\mu\) when solving for \(\vec{A}\).

Note

Application examples:

• cvp-i/, mms-v/, mms-vt-i/, ssol-i/, ssol-ii/, ssol-iii/.

Definition at line 318 of file static_vector_solver_i.hpp.

Constructor & Destructor Documentation

Solver1()

template<int dim, int stage = 1>

StaticVectorSolver::Solver1< dim, stage >::Solver1 ( unsigned int  p, unsigned int  mapping_degree, unsigned int  type_of_pde_rhs = 3, double  eta_squared = 0.0, std::string  fname = "data", const Function< dim > *  exact_solution = nullptr, bool  print_time_tables = false, bool  project_exact_solution = false, bool  write_higher_order_cells = false  )

inline

The only constructor.

Parameters

[in]	p	- Degree of the FE_Nedelec finite elements.
[in]	mapping_degree	- The degree of the interpolating polynomials used for mapping. Setting it to 1 will do in the most of the cases. Note, that it makes sense to attach a meaningful manifold to the triangulation if this parameter is greater than 1.
[in]	type_of_pde_rhs	- Defines how this class interprets the data provided by StaticVectorSolver::PdeRhs, see above.
[in]	eta_squared	- The gauging parameter, \(\eta^2\), in the partial differential equation.
[in]	fname	- The name of the output files without extension. Names of the output files will be generated by appending simulation conditions to this string.
[in]	exact_solution	- Points to an object that describes the exact solution to the problem. It is needed for calculating error norms. It is a responsibility of the user to make sure that the object exists at the time of the execution of run() or compute_error_norms().
[in]	print_time_tables	- If true, prints time tables on the screen.
[in]	project_exact_solution	- If true, projects the exact solution onto the space spanned by the Nedelec finite elements and saves the result into the output file next to the solution. This may be useful for debugging purposes as a comparison between the projected exact solution and the solution to the boundary value problem can yield a hint on where to search for bugs.
[in]	write_higher_order_cells	- Switches between the two modes of operation of the save() function, see the description of save().

Definition at line 352 of file static_vector_solver_i.hpp.

Member Function Documentation

fill_dirichlet_stack()

template<int dim, int stage = 1>

virtual void StaticVectorSolver::Solver1< dim, stage >::fill_dirichlet_stack ( )

pure virtual

Initializes the data member StaticVectorSolver::Solver1::dirichlet_stack.

This function must be overridden by the user. It must initialize the stack of the Dirichlet boundary conditions. For example:

using namespace dealii;
 
const types::boundary_id boundary_id_1 = 1;
const types::boundary_id boundary_id_2 = 3;
 
const DirichletFunction1<dim> dirichlet_bc_1;
const DirichletFunction2<dim> dirichlet_bc_2;
 
template<>
void SolverMyProblem<dim>::fill_dirichlet_stack()
{
  Solver<dim, stage>::dirichlet_stack =
    {{boundary_id_1, & dirichlet_bc_1},
     {boundary_id_2, & dirichlet_bc_2}};
}

The boundary IDs must be odd numbers, see above the convention on the boundary IDs.

make_mesh()

template<int dim, int stage = 1>

virtual void StaticVectorSolver::Solver1< dim, stage >::make_mesh ( )

pure virtual

Initializes the data member StaticVectorSolver::Solver1::triangulation.

This function must be overridden by the user. It must generate the mesh, label the boundaries, and, if necessary, assign user IDs. The mesh must be stored in the data member of this class StaticVectorSolver::Solver1::triangulation.

project_exact_solution_fcn()

template<int dim, int stage>

void StaticVectorSolver::Solver1< dim, stage >::project_exact_solution_fcn

Projects exact solution.

The mesh and the finite elements, are the same as are used for the numerical solution of the boundary value problem. The exact solution will be saved in the output file next to the numerical solution to the boundary value problem. This function works properly only if the exact solution is submitted to the constructor via the input parameter exact_solution and project_exact_solution=true.

Definition at line 1245 of file static_vector_solver_i.hpp.

run()

template<int dim, int stage = 1>

void StaticVectorSolver::Solver1< dim, stage >::run ( )

inline

Runs the simulation.

Executes the following member functions in a proper order: make_mesh(), fill_dirichlet_stack(), setup(), assemble(), solve(), project_exact_solution_fcn(), compute_error_norms(), save().

Definition at line 576 of file static_vector_solver_i.hpp.

save()

template<int dim, int stage>

void StaticVectorSolver::Solver1< dim, stage >::save

Saves simulation results into a vtk or vtu file.

The following data are saved:

• The calculated potential under the name "VectorField".
• The \(L^2\) error norm associated with the calculated potential under the name "L2norm". One value per mesh cell is saved.
• The \(L^{\infty}\) error norm associated with the calculated potential under the name "LinftyNorm". One value per mesh cell is saved.
• The exact solution expressed as a linear combination of the shape functions of the FE_Nedelec finite elements is saved under the name "VectorFieldExact". The "VectorField" and "VectorFieldExact" are modeled by exactly the same finite elements.

The "L2norm", "LinftyNorm", and "VectorFieldExact" are saved only if an exact solution is submitted to the constructor. Moreover, "VectorFieldExact" is calculated and saved only if project_exact_solution=true.

If write_higher_order_cells = false, the name of the file is computed by appending ".vtk" to the string contained by the parameter fname passed to the constructor. The file can be inspected with a help of VisIt or Paraview. Higher-order cells are saved as regular quadrilaterals and hexahedra. If write_higher_order_cells = true, the name of the file is computed by appending ".vtu" to the string contained by the parameter fname. The data is saved preserving the higher-order cells. The file can be viewed with a help of Paraview version 5.5.0 or higher. VisIt cannot load higher-order cells.

Definition at line 1262 of file static_vector_solver_i.hpp.

save_matrix_and_rhs_to_csv()

template<int dim, int stage>

void StaticVectorSolver::Solver1< dim, stage >::save_matrix_and_rhs_to_csv

(

std::string

fname

)

const

Saves the system matrix and the right-hand side into a csv file.

All the zeros are included into the csv files. This is a very dumb and inefficient way of saving sparse matrices. On the positive side - it is very easy and straightforward to read the csv files. This function may be useful for debugging. One can assemble the system on a coarse mesh (so there are a few mesh cells and the system matrix is small) and export the system matrix together with the right-hand side into another program such as Matlab or GNU Octave for an analysis.

Parameters

[in]

fname

- A stem of the names of the output files. The matrix will be saved into fname_matrix.csv file. The right-hand side will be save into fname_rhs.csv file.

Definition at line 1315 of file static_vector_solver_i.hpp.

setup()

template<int dim, int stage>

void StaticVectorSolver::Solver1< dim, stage >::setup

Initializes system matrix and the right-hand side vector.

Initialises StaticVectorSolver::Solver1::system_matrix, StaticVectorSolver::Solver1::system_rhs and other arrays. Applies the Dirichlet boundary conditions by constraining the system matrix. Distributes degrees of freedom.

Definition at line 765 of file static_vector_solver_i.hpp.

---

The documentation for this class was generated from the following file:

• static/shared/include/static_vector_solver_i.hpp