1518 *
const unsigned int mapping_degree = 2;
1519 *
const unsigned int fe_degree = 0;
1521 *
const double eta_squared = 0.0;
1523 *
const unsigned int n_threads_max = 0;
1525 *
const bool project_exact_solution =
false;
1531 * <a name=
"step_98-Convergencetable"></a>
1532 * <h3>Convergence table</h3>
1533 * The following
class describes a convergence table. The convergence tables are
1534 * saved on disk in TeX format.
1540 *
MainOutputTable() =
delete;
1542 *
MainOutputTable(
const unsigned int dim)
1547 *
void save(
const std::string &file_name)
1566 *
std::ofstream ofs(file_name +
".tex");
1571 *
const unsigned int dim;
1577 * <a name=
"step_98-Equations"></a>
1578 * <h3>Equations</h3>
1579 * The following
namespace contains closed-form analytical expressions
1580 *
for @f$T@f$, @f$\vec{
J}_f@f$, @f$\vec{
A}@f$, @f$B@f$, mentioned in the introduction to
this
1584 *
namespace ExactSolutions
1589 * The following function describes the free-current density, @f$\vec{
J}_f@f$,
1590 *
inside the current region. The current density in
this tutorial is
1591 * implemented by `ExactSolutions::FreeCurrentDensity`
class and by
1592 * `SolverT::Solver::free_current_density` member function . There is a
1593 * subtle difference in how these two implementation compute the free-current
1594 * density. Both classes, however, utilize the same expression
for the
1595 * free-current density. This function describes the expression.
1606 * The following
class implements the closed-form analytical
1607 * [expression](@ref Step98_Equation_1_Jf)
for the free-current density,
1608 * @f$\vec{
J}_f@f$, in the entire domain. The free-current density is computed
1609 * purely on the basis of the spatial coordinates of the field
point. In go
1610 * coordinates, out comes the free-current density. The information on the
1611 * material ID of the mesh cells and any other information on the mesh is
1612 * ignored. This function is used
for computing @f$L_2@f$ error norms and
for
1613 * computing the projected exact solution. The @f$\vec{
J}_f@f$ on the right-hand
1614 * side of the div-grad equation is implemented by the member function
1615 * `SolverT::Solver::free_current_density`
1618 *
class FreeCurrentDensity :
public Function<2>
1621 *
FreeCurrentDensity()
1630 *
ExcDimensionMismatch(
values.size(), p.size()));
1632 *
for (
unsigned int i = 0; i <
values.size(); i++)
1634 *
const double r = p[i].norm();
1636 *
if ((r >= Settings::a2) && (r <= Settings::b2))
1639 *
volume_free_current_density(p[i], Settings::K0);
1652 * The following
class implements the closed-form analytical
1653 * [expression](@ref Step98_Equation_4_A)
1654 *
for the magnetic vector potential, @f$\vec{
A}@f$.
1657 *
class MagneticVectorPotential :
public Function<2>
1660 *
MagneticVectorPotential()
1669 *
ExcDimensionMismatch(
values.size(), p.size()));
1671 *
for (
unsigned int i = 0; i <
values.size(); i++)
1673 *
const double r = p[i].norm();
1680 *
A = (mu_1 * K0 / 4.0) * (b2 * b2 - a2 * a2);
1684 *
A = (mu_0 * K0 / 4.0) * (b2 * b2 - a2 * a2) *
1685 *
(r + b1 * b1 * (mu_r - 1.0) / r) / r;
1689 *
A = (mu_0 * K0 / 2.0) *
1690 *
(b2 * b2 * r / 2.0 -
std::pow(r, 3) / 4.0 +
1691 *
(a2 * (b2 * b2 - a2 * a2) *
1692 *
(a2 + b1 * b1 * (mu_r - 1.0) / a2) / 2.0 -
1699 *
A = (mu_0 * K0 / 2.0) * b2 *
1701 *
(a2 * (b2 * b2 - a2 * a2) *
1702 *
(a2 + b1 * b1 * (mu_r - 1.0) / a2) / 2.0 -
1708 *
values[i][0] = -
A * p[i][1];
1716 * The following
class implements the closed-form analytical
1717 * [expression](@ref Step98_Equation_3_T)
for
1718 * the current vector potential, @f$T@f$.
1721 *
class CurrentVectorPotential : public
Function<2>
1724 *
CurrentVectorPotential()
1737 *
ExcMessage(
"This line is to avoid compiler warnings."));
1739 *
for (
unsigned int i = 0; i <
values.size(); i++)
1741 *
const double r = p[i].norm();
1746 *
values[i] = K0 * (b2 * b2 - a2 * a2) / 2.0;
1750 *
values[i] = -K0 * (r * r - b2 * b2) / 2.0;
1762 * The following
class implements the closed-form analytical
1763 * [expression](@ref Step98_Equation_2_B)
1764 *
for the magnetic field, @f$B@f$.
1767 *
class MagneticField : public
Function<2>
1783 *
ExcMessage(
"This line is to avoid compiler warnings."));
1785 *
for (
unsigned int i = 0; i <
values.size(); i++)
1787 *
const double r = p[i].norm();
1792 *
values[i] = mu_1 * K0 * (b2 * b2 - a2 * a2) / 2.0;
1796 *
values[i] = mu_0 * K0 * (b2 * b2 - a2 * a2) / 2.0;
1800 *
values[i] = mu_0 * K0 * (b2 * b2 - r * r) / 2.0;
1815 * <a name=
"step_98-BaseSolver"></a>
1816 * <h3>Base Solver</h3>
1820 * As discussed above, the following
namespace aggregates the code common to all
1821 * four solvers used in the tutorial. All four solvers are derived from the
1822 * `BaseSolver`
class.
1825 *
namespace BaseClasses
1829 * The computation of the integrals of the functionals is delegated to the
1830 * derived classes. The function that computes the integrals,
1831 * `BaseSolver::system_matrix_local`, is
virtual and must be overridden by
1832 * the derived classes. However, the objects of the
types FEValues are
1833 * initialized at the
level of the `BaseSolver`
class together with
1834 * `AssemblyScratchData`. For the initialization to work properly, the
1835 * `BaseSolver`
class must know which cell
data to compute
for a particular
1836 * implementation of the solver down the hierarchy. This information is
1837 * communicated to the `BaseSolver` by passing an argument of the type
1838 * `UpdateFlagsCollection` to the constructor. In all solvers, with exception
1839 * of `SolverT`, we have two
types of finite elements. One type of finite
1840 * elements models the solution to the partial differential equation. Another
1841 * models the physical quantity on the right-hand side of the partial
1842 * differential equation. The `solution_update_flags`
data member below
1843 * contains the flags
for updating the
values of the finite elements that
1844 * model the solution. The `rhs_update_flags`
data member contains the flags
1845 *
for updating the
values of the finite elements that model the physical
1846 * quantity on the right-hand side of the partial differential equation.
1849 *
struct UpdateFlagsCollection
1857 * Each iteration of the program consists of
1858 * [four stages](@ref Step98_FourStages). Each stage utilizes one solver. All
1859 * four solvers used in
this tutorial are derived from the following
class.
1860 * The solver used in the
first stage loads the mesh and refines it
if
1861 * necessary. The solvers in the other three stages reuse the mesh prepared at
1862 * the
first stage. Furthermore, the solver in the
first stage expects a
1863 * closed-form analytical expression on the right-hand side of the partial
1864 * differential equation. Each solver in the other three stages expects a
1865 * potential computed at one of the preceding stages, i.e., a field in a form
1866 * of linear superposition of the shape
functions. At the top of the following
1867 *
class, we declare two constructors. The
first constructor must be used
for
1868 * constructing solver
for the
first stage. The
second constructor must be
1869 * used
for constructing the solvers
for the
second, third, and fourth stages.
1870 * The
second constructor has three extra arguments, `triangulation_rhs`,
1871 * `dof_handler_rhs`, and `solution_rhs`
for accommodating the mesh and the
1872 * potential computed at one of the preceding stages.
1876 * Following the constructors are the eight functions that implement various
1877 * steps typical
for every solver. The function `run` aggregates these steps.
1878 * This arrangement of functions is quite standard in deal.II, see @ref step_3
"step-3",
1879 *
for instance. Normally, the `setup` function begins by distributing the
1880 * dofs. The problem is: the type of the finite elements is not known in the
1881 * `BaseSolver` class. It is specified in the derived classes. Therefore, it
1882 * could be reasonable to make the setup function virtual as well. Instead,
1883 * we move the dof distribution code to the
end of the `make_mesh` function
1884 * which is, in fact, virtual and must be overridden in the derived classes
1889 * After that, we declare six get
functions that provide access to protected
1891 * `MagneticProblem::
run` function. The
data provided by these get-
functions
1892 * is used to fill the convergence tables and to pass the references to the
1893 * triangulation, the dof handler, and the dofs between solvers.
1899 *
BaseSolver(
const unsigned
int mapping_degree,
1900 *
const UpdateFlagsCollection update_flags_collection,
1907 *
const unsigned int stage,
1908 *
const unsigned int mapping_degree,
1909 *
const UpdateFlagsCollection update_flags_collection,
1910 *
const std::string &file_name,
1913 *
virtual void make_mesh() = 0;
1920 *
void output_results()
const;
1921 *
void compute_error_norms();
1922 *
void project_exact_solution_fcn();
1927 *
double get_L2_norm()
const;
1928 *
unsigned int get_n_cells()
const;
1936 * We
begin the `
protected` section of the `BaseSolver`
class by declaring
1937 * three
data members that store the input from
one of the preceding
1938 * solvers. If there is no preceding solver and these
data is not provided,
1939 * i.e.,
first of the two constructors above has been used, these three
data
1940 * members
point to `triangulation`, `dof_handler`, and `solution`, of the
1945 * Following are the three
data members that store the triangulation, dof
1946 * handler, and dofs vector of the current solver. In the
case of the
1947 *
first-stage solver, the `triangulation`
data member stores the loaded and
1948 * refined mesh. This
data member is not used in the
case of the solvers of
1949 * the
second, third, and fourth stages. The `dof_handler` and `solution`
1950 * represent the result of the solver, i.e., the computed potential or
1955 * Next, we declare four
data members that describe the system of linear
1956 * equations to be solved by the linear solver. The components of the system
1957 *
matrix and that of the right-hand side are computed by the `
assemble`
1958 * function. The
affine constraints are used to
apply the Dirichlet boundary
1959 * conditions and to distribute the local (cell specific) system
matrix and
1960 * right-hand side to `system_matrix` and `system_rhs`. The are no hanging
1961 * nodes and hanging node constraints in
this program. The last
data member
1962 * in
this block describes the
dynamic sparsity pattern. The tutorial @ref step_2
"step-2"
1963 * discusses the rationale behind the
dynamic sparsity pattern.
1967 * The next block contains two
data members related to the exact solution.
1968 * The
data member `exact_solution` points to the closed-form analytical
1969 * solution the solver attempts to compute.
1970 * If `Settings::project_exact_solution=
true`, the exact solution is
1971 * projected onto a proper function space and the
data member
1972 * `projected_exact_solution` is populated by the dofs of the projected
1973 * exact solution. The corresponding dof handler is `dof_handler`. Together
1974 * `dof_handler` and `projected_exact_solution` constitute the field
1975 * function which describes the exact solution. It is saved into the `.vtu`
1976 * file next to the solution and the @f$L_2@f$ error
norm.
1980 * The next four
data members simply store the
data supplied as arguments
1981 * to the constructor. The `stage`
data member stores the number of the
1982 * current stage. The `mapping_degree`
data member contains the degree of
1983 * mapping from the
reference cell to a mesh cell and back. The
1984 * `update_flags_collection` contains the information on which finite
1985 * element
values must be computed
for each cell. The names of the output
1986 * files are derived by appending strings to `file_name`.
1990 * The next block contains two
data members computed by the function
1991 * `BaseSolver::compute_error_norms`. The `L2_per_cell`
data member
1992 * contains
one value of the @f$L^2@f$ error
norm per mesh cell. It is saved
1993 * into the `.vtu` file next to the solution. The `
L2_norm`
data member
1994 * contains
one value of the @f$L^2@f$ error
norm per mesh. It is reported in
1995 * the convergence table.
2015 *
const unsigned int stage;
2016 *
const unsigned int mapping_degree;
2017 *
const UpdateFlagsCollection update_flags_collection;
2018 *
const std::string file_name;
2025 * The program utilizes the
WorkStream technology. The @ref step_9
"step-9" tutorial
2026 * does a much better job of explaining the workings of
WorkStream.
2027 * Reading the @ref workstream_paper
"WorkStream paper" is recommended.
2028 * In very simple terms, the workings of the
WorkStream can be envisioned as
2029 * the following. Let us assume we have a task of computing components of
2030 * the system
matrix, @f$A_{ij}@f$, and the right-hand side vector, @f$b_i@f$.
2031 * Simply put, we need to fill in the
matrix `system_matrix` and vector
2032 * `system_rhs`. This is a big task as the number of dofs is large. The
2033 * idea is to
split the big task on a number of small tasks and feed them
2034 * to multiple threads to speed up the calculation process. Each small
2035 * task consists of computing the contributions of a single mesh cell to
2036 * `system_matrix` and `system_rhs`. These contributions are stored
2037 * temporary in `
cell_matrix` and `cell_rhs`
for each cell. These
2038 * contributions are then copied to `system_matrix` and `system_rhs`.
2039 *
WorkStream creates and schedules the small tasks and takes care
2040 * of copying `
cell_matrix` and `cell_rhs` to `system_matrix` and
2041 * `system_rhs`. The rest of the declarations in the `
protected` section
2042 * of the `BaseSolver`
class help to communicate to
WorkStream information
2043 * which is necessary
for its operation.
2047 * First, the `CellIteratorPair` type is declared in the block of code
2048 * below. The solvers in the
second, third, and fourth stages use two dof
2049 * handlers, `dof_handler` and `dof_handler_rhs`. The
WorkStream needs to
2050 * walk through the two dof handlers synchronously. For
this purpose we
2051 * pair two active cell iterators (one from `dof_handler`, another from
2052 * `dof_handler_rhs`). For that we need the `CellIteratorPair` type. The
2053 * solver at the
first stage (a solver constructed by invoking the
first
2054 * constructor above) uses only
one dof handler. In
this case the
2055 * constructor makes `dof_handler_rhs` to
reference `dof_handler`. In
2056 * effect, both iterators of the tuple will iterate the same dof handler.
2057 * In the
case of the
first-stage solver we will use only the
first
2062 * Next, we declare the `AssemblyScratchData` type. An
object of
this type
2063 * contains the relevant information on the current mesh cell which is used
2064 * as an input
for computing the components of the system
matrix and the
2065 * right-hand side.
WorkStream creates an
object of
this type and passes
2066 * it to the function `system_matrix_local` which, in turn, computes
2071 * Next, the type `AssemblyCopyData` is declared. An objects of
this type
2072 * contains the cell specific contributions to the system
matrix and the
2073 * right-hand side. The
WorkStream creates
object of
this type and passes
2074 * it to function `system_matrix_local` along with an
object of the type
2075 * `AssemblyScratchData`. The function `system_matrix_local`, in turn,
2076 * takes input
data form the
object of the type `AssemblyScratchData`,
2077 * computes the relevant integrals and places the result into the
object
2078 * of the type `AssemblyCopyData`. The
WorkStream copies the content of
2079 * the `AssemblyCopyData` object, `
cell_matrix` and `cell_rhs`, into
2080 * `system_matrix` and `system_rhs`. The
data member
2081 * `AssemblyCopyData::local_dof_indices` contains the indices of global
2082 * components to which cell-specific
data must be copied.
2086 *
WorkStream calls the function `system_matrix_local` to compute the
2087 * cell-specific components. Likewise,
WorkStream uses calls to
2088 * `copy_local_to_global` to
copy the cell-specific
data into the system
2090 * `copy_local_to_global` are declared last.
2093 *
using IteratorTuple =
2094 *
std::tuple<typename DoFHandler<2>::active_cell_iterator,
2099 *
struct AssemblyScratchData
2104 *
const unsigned int mapping_degree,
2105 *
const UpdateFlagsCollection update_flags_collection,
2106 *
const unsigned int stage);
2108 *
AssemblyScratchData(
const AssemblyScratchData &scratch_data);
2115 *
const unsigned int dofs_per_cell;
2116 *
const unsigned int n_q_points;
2118 *
std::vector<double> permeability_list;
2119 *
std::vector<double> values_list_rhs;
2120 *
std::vector<Tensor<1, 2>> vectors_list_rhs;
2121 *
std::vector<std::vector<Tensor<1, 2>>> vectors_vectors_list_rhs;
2127 *
struct AssemblyCopyData
2131 *
std::vector<types::global_dof_index> local_dof_indices;
2134 *
virtual void system_matrix_local(
const CellIteratorPair &IP,
2135 *
AssemblyScratchData &scratch_data,
2136 *
AssemblyCopyData ©_data) = 0;
2138 *
void copy_local_to_global(
const AssemblyCopyData ©_data);
2143 * The following are the implementations of the two constructors
2144 * of the `BaseSolver`
class.
2147 *
BaseSolver::BaseSolver(
const unsigned int mapping_degree,
2148 *
const UpdateFlagsCollection update_flags_collection,
2149 *
const std::string &file_name,
2151 *
: triangulation_rhs(triangulation)
2152 *
, dof_handler_rhs(dof_handler)
2153 *
, solution_rhs(solution)
2154 *
, exact_solution(exact_solution)
2156 *
, mapping_degree(mapping_degree)
2157 *
, update_flags_collection(update_flags_collection)
2158 *
, file_name(file_name)
2164 *
const unsigned int stage,
2165 *
const unsigned int mapping_degree,
2166 *
const UpdateFlagsCollection update_flags_collection,
2167 *
const std::string &file_name,
2169 *
: triangulation_rhs(triangulation_rhs)
2170 *
, dof_handler_rhs(dof_handler_rhs)
2171 *
, solution_rhs(solution_rhs)
2172 *
, exact_solution(exact_solution)
2174 *
, mapping_degree(mapping_degree)
2175 *
, update_flags_collection(update_flags_collection)
2176 *
, file_name(file_name)
2181 * The following function applies the Dirichlet boundary condition, sets
2182 * up a sparsity pattern, and initializes the vectors and matrices. It is
2183 * common
for the setup function to distribute the dofs. The type of the
2184 * finite elements, however, is not known at the
level of `BaseSolver`
class.
2185 * The type of the finite elements is chosen in the derived classes. For
2186 *
this reason, the task of distributing the dofs is shifted to the
end of
2187 * the `make_mesh` function which is a
virtual function.
2191 * The Dirichlet boundary condition must be enforced only in the
first
2192 * [stage](@ref Step98_FourStages)
2193 * where the div-grad equation is solved
for the current vector potential,
T.
2194 * For
this reason we have `
if (stage == 1)` filter in the beginning of the
2195 * function. The boundary
value problem
for the curl-curl equation utilizes
2196 * the Neumann boundary condition. It is a natural boundary condition. It is
2197 * enforced by minimization of the functional. The two projectors,
2198 * @f$T \rightarrow \vec{
J}_f@f$ and @f$\vec{
A} \rightarrow B@f$, use no boundary
2202 *
void BaseSolver::setup()
2204 *
constraints.clear();
2211 *
Settings::outer_boundary_id,
2215 *
constraints.close();
2219 *
sparsity_pattern.copy_from(dsp);
2221 *
system_matrix.reinit(sparsity_pattern);
2222 *
solution.reinit(dof_handler.n_dofs());
2223 *
system_rhs.reinit(dof_handler.n_dofs());
2225 *
if (Settings::project_exact_solution && exact_solution)
2226 *
projected_exact_solution.reinit(dof_handler.n_dofs());
2228 *
if (exact_solution)
2229 *
L2_per_cell.reinit(triangulation.n_active_cells());
2234 * Formally, the following function assembles the system of linear equations.
2235 * In reality, however, it just spells all the magic words to get the
2236 *
WorkStream going. The interesting part, i.e., computing the components of
2237 * the system
matrix and the right-hand side, happens in the function
2238 * `system_matrix_local` which is a
virtual function. That is to say, the
2239 * recipe
for the functional is implemented in the derived
class by overriding
2240 * the
virtual function `system_matrix_local`.
2243 *
void BaseSolver::assemble()
2246 *
dof_handler_rhs.begin_active()}),
2248 *
{dof_handler.end(), dof_handler_rhs.end()}),
2250 *
&BaseSolver::system_matrix_local,
2251 *
&BaseSolver::copy_local_to_global,
2252 *
AssemblyScratchData(dof_handler,
2256 *
update_flags_collection,
2258 *
AssemblyCopyData());
2263 * The following are the implementation of the constructors of the
2264 * `AssemblyScratchData`. The
first constructor initializes the scratch
2265 *
data from the input parameters. The
second - from another
object of
2266 * the same type, i.e., a
copy constructor.
2269 *
BaseSolver::AssemblyScratchData::AssemblyScratchData(
2273 *
const unsigned int mapping_degree,
2274 *
const UpdateFlagsCollection update_flags_collection,
2275 *
const unsigned int stage)
2276 *
: mapping(mapping_degree)
2277 *
, fe_values(mapping,
2278 *
dof_handler.get_fe(),
2279 *
QGauss<2>((stage == 1) ? (dof_handler.get_fe().degree + 1) :
2280 *
(dof_handler.get_fe().degree + 2)),
2281 *
update_flags_collection.solution_update_flags)
2282 *
, fe_values_rhs(mapping,
2283 *
dof_handler_rhs.get_fe(),
2284 *
QGauss<2>((stage == 1) ? (dof_handler.get_fe().degree + 1) :
2285 *
(dof_handler.get_fe().degree + 2)),
2286 *
update_flags_collection.rhs_update_flags)
2287 *
, dofs_per_cell(fe_values.dofs_per_cell)
2288 *
, n_q_points(fe_values.get_quadrature().
size())
2289 *
, permeability_list(n_q_points)
2290 *
, values_list_rhs(n_q_points)
2291 *
, vectors_list_rhs(n_q_points)
2292 *
, vectors_vectors_list_rhs(n_q_points,
std::vector<
Tensor<1, 2>>(2))
2293 *
, dof_handler_rhs(dof_handler_rhs)
2294 *
, dofs_rhs(dofs_rhs)
2297 *
BaseSolver::AssemblyScratchData::AssemblyScratchData(
2298 *
const AssemblyScratchData &scratch_data)
2299 *
: mapping(scratch_data.mapping.
get_degree())
2300 *
, fe_values(mapping,
2301 *
scratch_data.fe_values.get_fe(),
2302 *
scratch_data.fe_values.get_quadrature(),
2304 *
, fe_values_rhs(mapping,
2305 *
scratch_data.fe_values_rhs.get_fe(),
2306 *
scratch_data.fe_values_rhs.get_quadrature(),
2308 *
, dofs_per_cell(fe_values.dofs_per_cell)
2309 *
, n_q_points(fe_values.get_quadrature().
size())
2310 *
, permeability_list(n_q_points)
2311 *
, values_list_rhs(n_q_points)
2312 *
, vectors_list_rhs(n_q_points)
2313 *
, vectors_vectors_list_rhs(n_q_points,
std::vector<
Tensor<1, 2>>(2))
2314 *
, dof_handler_rhs(scratch_data.dof_handler_rhs)
2315 *
, dofs_rhs(scratch_data.dofs_rhs)
2320 * The following function copies the components of a cell
matrix and a cell
2321 * right-hand side into the system
matrix, @f$A_{ij}@f$, and the system right-hand
2325 *
void BaseSolver::copy_local_to_global(
const AssemblyCopyData ©_data)
2327 *
constraints.distribute_local_to_global(copy_data.
cell_matrix,
2328 *
copy_data.cell_rhs,
2329 *
copy_data.local_dof_indices,
2336 * The following function solves the system of linear equations.
2337 * The stopping condition
for the iteration algorithm is
2338 * @f$\|\boldsymbol{
b}-\boldsymbol{
A}\boldsymbol{c}\|<10^{-8}\|\boldsymbol{
b}\|@f$.
2339 * The maximum number of iteration steps is
set to `system_rhs.size()` as
2340 * the conjugate
gradient algorithm is supposed to find the solution in at
2341 * most @f$m@f$ steps
for an @f$m \times m@f$ system
matrix. This function also
2342 * distributes constraints. The constraints are only used to enforce the
2343 * Dirichlet boundary condition.
2346 *
void BaseSolver::solve()
2348 *
SolverControl control(system_rhs.size(), 1e-8 * system_rhs.l2_norm());
2354 *
preconditioner.
initialize(system_matrix, 1.2);
2356 *
cg.solve(system_matrix, solution, system_rhs, preconditioner);
2358 *
constraints.distribute(solution);
2363 * The following two
functions compute the error norms and
project the exact
2367 *
void BaseSolver::compute_error_norms()
2369 *
if (exact_solution)
2377 *
dof_handler.get_fe(0).degree + 4),
2386 *
void BaseSolver::project_exact_solution_fcn()
2388 *
if (Settings::project_exact_solution && exact_solution)
2392 *
constraints_empty.
clear();
2396 *
constraints_empty.close();
2400 *
constraints_empty,
2402 *
(dof_handler.get_fe().degree + 1) :
2403 *
(dof_handler.get_fe().degree + 2)),
2405 *
projected_exact_solution);
2411 * The following function saves the solution and the @f$L_2@f$ error
norm into a
2412 * `.vtu` file. If `Settings::project_exact_solution =
true`, the projected
2413 * exact solution is saved as well.
2416 *
void BaseSolver::output_results() const
2418 *
const bool fe_is_vector = ((stage == 2) || (stage == 3));
2420 *
const std::string name = (fe_is_vector ?
"VectorField" :
"ScalarField");
2421 *
const unsigned int components = (fe_is_vector ? 2 : 1);
2423 *
component_interpretation =
2428 *
std::vector<std::string> solution_names(components, name);
2429 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
2430 *
data_component_interpretation(components, component_interpretation);
2437 *
data_component_interpretation);
2439 *
if (exact_solution)
2441 *
data_out.add_data_vector(L2_per_cell,
"L2norm");
2443 *
if (Settings::project_exact_solution)
2445 *
std::vector<std::string> solution_names_exact(components,
2447 *
data_out.add_data_vector(dof_handler,
2448 *
projected_exact_solution,
2449 *
solution_names_exact,
2450 *
data_component_interpretation);
2456 *
data_out.set_flags(flags);
2460 *
data_out.build_patches(mapping,
2461 *
dof_handler.get_fe(0).degree + 2,
2464 *
std::ofstream ofs(file_name +
".vtu");
2465 *
data_out.write_vtu(ofs);
2470 * The following function clears the memory
for the next solver.
2473 *
void BaseSolver::clear()
2475 *
system_matrix.clear();
2476 *
system_rhs.reinit(0);
2485 *
void BaseSolver::run()
2491 *
compute_error_norms();
2492 *
project_exact_solution_fcn();
2499 * The following is the straightforward implementation of the get
functions.
2502 *
double BaseSolver::get_L2_norm() const
2507 *
unsigned int BaseSolver::get_n_cells() const
2509 *
return triangulation_rhs.n_active_cells();
2514 *
return dof_handler.n_dofs();
2519 *
return triangulation_rhs;
2524 *
return dof_handler;
2536 * <a name=
"step_98-SolverT"></a>
2537 * <h3>Solver -
T</h3>
2541 * The following
namespace contains the code related to the computation of the
2542 * current vector potential, @f$T@f$.
2547 *
using namespace BaseClasses;
2551 * We derive the solver from the `BaseSolver`
class. What is left to
do is
2552 * to initialize the `BaseSolver`,
override the two
virtual functions
2553 * (`make_mesh` and `system_matrix_local`), and implement the function
2554 * `free_current_density`. The function `free_current_density` implements
2555 * the closed-form analytical expression
for @f$\vec{
J}_f@f$ on the right-hand
2556 * side of the [div-grad equation](@ref Step98_FourStages).
2559 *
class Solver : public BaseSolver
2562 *
Solver() = delete;
2563 *
Solver(
const unsigned int p,
2564 *
const unsigned int r,
2565 *
const unsigned int mapping_degree,
2566 *
const std::string &file_name =
"data",
2570 *
virtual void make_mesh() override final;
2572 *
system_matrix_local(const CellIteratorPair &IP,
2573 *
AssemblyScratchData &scratch_data,
2574 *
AssemblyCopyData ©_data) override final;
2577 *
const
unsigned int refinement_parameter;
2579 *
void free_current_density(const
std::vector<
Point<2>> &p,
2580 *
const
types::material_id material_id,
2581 *
std::vector<
Tensor<1, 2>> &values) const;
2586 * Following is the implementation of the constructor.
2587 * We use the
first constructor of the `BaseSolver` class as
2588 * the solver is used at the
2589 * [
first stage](@ref Step98_FourStages). By looking at the
2590 * [expressions](@ref Step98_Numerical_Recipe_T)
2591 *
for @f$A_{ij}@f$ and @f$b_i@f$ we can conclude that to compute them we need
2593 * weights multiplied by the Jacobian
determinant(`JxW`). The quadrature
2594 * points are needed to sample the closed-form analytical expression
for
2595 * @f$\vec{
J}_f@f$. Accordingly, we use the update flags
2597 *
for the
FE_Q finite elements. This time we
do not use the finite elements
2598 * that model the potential on the right-hand side of the equation, so we
2602 *
Solver::Solver(
const unsigned int p,
2603 *
const unsigned int r,
2604 *
const unsigned int mapping_degree,
2605 *
const std::string &file_name,
2607 *
: BaseSolver(mapping_degree,
2615 *
, refinement_parameter(r)
2620 * The following function loads the mesh, creates manifolds, bounds the
2621 * manifolds to the manifold IDs, refines the mesh, and distributes the dofs.
2622 * Note that we are allowed to create the manifolds locally in the function as
2623 * the triangulation
object keeps copies of the manifolds, see @ref step_65
"step-65". Also
2624 * recall that we have shifted the task of distributing the dofs from the
2625 * `setup` function to the `make_mesh` function to evade the necessity to make
2626 * the `setup` function
virtual. This allows us to keep
one `setup` function
2627 * in the `BaseSolver`
class that serves the needs of all derived classes.
2630 *
void Solver::make_mesh()
2635 *
gridin.read_msh(
"circle.msh");
2643 *
transfinite_manifold.
initialize(triangulation);
2644 *
triangulation.set_manifold(transfinite_interpolation_manifold_id,
2645 *
transfinite_manifold);
2647 *
triangulation.refine_global(refinement_parameter);
2649 *
dof_handler.reinit(triangulation);
2650 *
dof_handler.distribute_dofs(fe);
2655 * The following function assembles a fraction of
2656 * [the system
matrix and the system right-hand side](@ref Step98_Numerical_Recipe_T)
2657 * related to a single cell. These fractions are
2658 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied to
2659 * `system_matrix` and `system_rhs` by
WorkStream.
2662 *
void Solver::system_matrix_local(
const CellIteratorPair &IP,
2663 *
AssemblyScratchData &scratch_data,
2664 *
AssemblyCopyData ©_data)
2668 *
copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
2669 *
scratch_data.dofs_per_cell);
2671 *
copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
2673 *
copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
2677 *
scratch_data.fe_values.reinit(cell);
2679 *
Solver::free_current_density(scratch_data.fe_values.get_quadrature_points(),
2680 *
cell->material_id(),
2681 *
scratch_data.vectors_list_rhs);
2683 *
for (
unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index)
2685 *
for (
unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i)
2687 *
for (
unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j)
2689 *
copy_data.cell_matrix(i, j) +=
2690 *
(scratch_data.fe_values[se].gradient(
2692 *
scratch_data.fe_values[se].gradient(
2695 *
scratch_data.fe_values.JxW(q_index);
2697 *
copy_data.cell_rhs(i) +=
2698 *
(scratch_data.vectors_list_rhs[q_index][0] *
2699 *
scratch_data.fe_values[se].gradient(i, q_index)[1] -
2700 *
scratch_data.vectors_list_rhs[q_index][1] *
2701 *
scratch_data.fe_values[se].gradient(i,
2703 *
scratch_data.fe_values.JxW(
2708 *
cell->get_dof_indices(copy_data.local_dof_indices);
2713 * The following function implements the closed form analytical
2714 * [expression](@ref Step98_Equation_1_Jf)
2715 *
for @f$\vec{
J}_f@f$ on the right-hand side of the
2716 * [div-grad equation](@ref Step98_PDE_T).
2719 *
void Solver::free_current_density(
const std::vector<
Point<2>> &p,
2724 *
ExcDimensionMismatch(p.size(),
values.size()));
2726 *
if ((material_id == Settings::material_id_free_space) ||
2727 *
(material_id == Settings::material_id_core))
2730 *
if (material_id == Settings::material_id_free_current)
2731 *
for (
unsigned int i = 0; i <
values.size(); i++)
2733 *
ExactSolutions::volume_free_current_density(p[i], Settings::K0);
2740 * <a name=
"step_98-ProjectorfromHgradtoHdiv"></a>
2741 * <h3>Projector from H(grad) to H(div)</h3>
2745 * The following
namespace contains all the code related to the computation of
2746 * the free-current density, @f$\vec{
J}_f@f$.
2749 *
namespace ProjectorHgradToHdiv
2751 *
using namespace BaseClasses;
2755 * We derive the solver from the `BaseSolver`
class. What is left to
do is
2756 * to initialize the `BaseSolver` and
override the two
virtual functions
2757 * (`make_mesh` and `system_matrix_local`).
2760 *
class Solver : public BaseSolver
2763 *
Solver() = delete;
2764 *
Solver(
const unsigned int p,
2766 *
const unsigned int mapping_degree,
2770 *
const std::string &file_name =
"data",
2774 *
virtual void make_mesh() override final;
2776 *
system_matrix_local(const CellIteratorPair &IP,
2777 *
AssemblyScratchData &scratch_data,
2778 *
AssemblyCopyData ©_data) override final;
2785 * Following is the implementation of the constructor.
2786 * We use the
second constructor of the `BaseSolver` class as
2787 * the solver is used at the
2788 * [
second stage](@ref Step98_FourStages). By looking at the
2789 * [expressions](@ref Step98_Numerical_Recipe_Jf)
2790 *
for @f$A_{ij}@f$ and @f$b_i@f$ we can conclude that to compute them we need
2791 *
values of the shape
functions and the quadrature weights multiplied by the
2795 * there is a numerically computed potential, @f$T@f$, on the right-hand side of
2796 * the [equation](@ref Step98_PDE_Jf). It is modeled by the
FE_Q finite elements.
2797 * To compute the right-hand side, we need gradients of the shape
functions.
2802 *
Solver::Solver(
const unsigned int p,
2803 *
const unsigned int mapping_degree,
2807 *
const std::string &file_name,
2809 *
: BaseSolver(triangulation_rhs,
2822 * At the
second stage we
do not
load the mesh. We reuse the mesh loaded
2823 * at the
first stage. Consequently, we just need to distribute the dofs.
2826 *
void Solver::make_mesh()
2828 *
dof_handler.reinit(triangulation_rhs);
2829 *
dof_handler.distribute_dofs(fe);
2834 * The following function assembles a fraction of
2835 * [the system
matrix and the system right-hand side](@ref Step98_Numerical_Recipe_Jf)
2836 * related to a single cell. These fractions are
2837 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied to
2838 * `system_matrix` and `system_rhs` by
WorkStream.
2841 *
void Solver::system_matrix_local(
const CellIteratorPair &IP,
2842 *
AssemblyScratchData &scratch_data,
2843 *
AssemblyCopyData ©_data)
2847 *
copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
2848 *
scratch_data.dofs_per_cell);
2850 *
copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
2852 *
copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
2858 *
scratch_data.fe_values.reinit(cell);
2859 *
scratch_data.fe_values_rhs.reinit(cell_rhs);
2861 *
scratch_data.fe_values_rhs.get_function_gradients(
2862 *
scratch_data.dofs_rhs, scratch_data.vectors_list_rhs);
2864 *
for (
unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index)
2866 *
for (
unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i)
2868 *
for (
unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j)
2870 *
copy_data.cell_matrix(i, j) +=
2871 *
scratch_data.fe_values[ve].value(i, q_index) *
2872 *
scratch_data.fe_values[ve].value(j,
2874 *
scratch_data.fe_values.JxW(q_index);
2877 *
copy_data.cell_rhs(i) +=
2878 *
(scratch_data.vectors_list_rhs[q_index][1] *
2879 *
scratch_data.fe_values[ve].value(i, q_index)[0] -
2880 *
scratch_data.vectors_list_rhs[q_index][0] *
2881 *
scratch_data.fe_values[ve].value(i, q_index)[1]) *
2882 *
scratch_data.fe_values.JxW(
2887 *
cell->get_dof_indices(copy_data.local_dof_indices);
2896 * <a name=
"step_98-SolverA"></a>
2897 * <h3>Solver -
A</h3>
2901 * The following
namespace contains all the code related to the computation of
2902 * the magnetic vector potential, @f$\vec{
A}@f$.
2907 *
using namespace BaseClasses;
2911 * We derive the solver from the `BaseSolver`
class. What is left to
do is
2912 * to initialize the `BaseSolver`,
override the two
virtual functions
2913 * (`make_mesh` and `system_matrix_local`), and implement the function
2917 *
class Solver : public BaseSolver
2920 *
Solver() = delete;
2921 *
Solver(
const unsigned int p,
2922 *
const unsigned int mapping_degree,
2926 *
const std::string &file_name =
"data",
2930 *
virtual void make_mesh() override final;
2932 *
system_matrix_local(const CellIteratorPair &IP,
2933 *
AssemblyScratchData &scratch_data,
2934 *
AssemblyCopyData ©_data) override final;
2938 *
void permeability(const
types::material_id material_id,
2939 *
std::vector<
double> &values) const;
2944 * Following is the implementation of the constructor.
2945 * We use the
second constructor of the `BaseSolver` class as
2946 * the solver is used at the
2947 * [third stage](@ref Step98_FourStages). By looking at the
2948 * [expressions](@ref Step98_Numerical_Recipe_A)
2949 *
for @f$A_{ij}@f$ and @f$b_i@f$ we can conclude that to compute them we need
2952 * finite elements. Accordingly, we use the update flags `
update_values`,
2954 * elements. This time there is a numerically computed potential, @f$T@f$, on
2955 * the right-hand side of the
2956 * [equation](@ref Step98_PDE_A). It is modeled by the
FE_Q finite elements.
2957 * To compute it, we need values of the shape
functions. Accordingly, we use
2961 *
Solver::Solver(
const unsigned int p,
2962 *
const unsigned int mapping_degree,
2966 *
const std::string &file_name,
2968 *
: BaseSolver(triangulation_rhs,
2982 * At the third stage we
do not
load the mesh. We reuse the mesh loaded at
2983 * the
first stage. Consequently, we just need to distribute the dofs.
2986 *
void Solver::make_mesh()
2988 *
dof_handler.reinit(triangulation_rhs);
2989 *
dof_handler.distribute_dofs(fe);
2994 * The following function assembles a fraction of
2995 * [the system
matrix and the system right-hand side](@ref Step98_Numerical_Recipe_A)
2996 * related to a single cell. These fractions are
2997 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied to
2998 * `system_matrix` and `system_rhs` by
WorkStream.
3001 *
void Solver::system_matrix_local(
const CellIteratorPair &IP,
3002 *
AssemblyScratchData &scratch_data,
3003 *
AssemblyCopyData ©_data)
3007 *
copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
3008 *
scratch_data.dofs_per_cell);
3010 *
copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
3012 *
copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
3018 *
scratch_data.fe_values.reinit(cell);
3019 *
scratch_data.fe_values_rhs.reinit(cell_rhs);
3021 *
Solver::permeability(cell->material_id(), scratch_data.permeability_list);
3023 *
scratch_data.fe_values_rhs.get_function_values(
3024 *
scratch_data.dofs_rhs, scratch_data.values_list_rhs);
3026 *
for (
unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index)
3028 *
for (
unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i)
3030 *
for (
unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j)
3032 *
copy_data.cell_matrix(i, j) +=
3033 *
(1.0 / scratch_data.permeability_list[q_index]) *
3034 *
(scratch_data.fe_values[ve].curl(
3036 *
scratch_data.fe_values[ve].curl(
3039 *
Settings::eta_squared *
3040 *
scratch_data.fe_values[ve].value(i,
3042 *
scratch_data.fe_values[ve].value(j, q_index)
3044 *
scratch_data.fe_values.JxW(q_index);
3046 *
copy_data.cell_rhs(i) +=
3047 *
(scratch_data.values_list_rhs[q_index] *
3048 *
scratch_data.fe_values[ve].curl(i, q_index)) *
3049 *
scratch_data.fe_values.JxW(q_index);
3053 *
cell->get_dof_indices(copy_data.local_dof_indices);
3058 * The following function implements the
3059 * [equation](@ref Step98_Equation_MU)
3064 *
std::vector<double> &values)
const
3066 *
if (material_id == Settings::material_id_core)
3067 *
std::fill(
values.begin(),
values.end(), Settings::mu_1);
3069 *
std::fill(
values.begin(),
values.end(), Settings::mu_0);
3076 * <a name=
"step_98-ProjectorfromHcurltoL2"></a>
3077 * <h3>Projector from H(curl) to
L2</h3>
3081 * The following
namespace contains all the code related to the computation of
3082 * the magnetic field, @f$B@f$.
3085 *
namespace ProjectorHcurlToL2
3087 *
using namespace BaseClasses;
3091 * We derive the solver from the `BaseSolver`
class. What is left to
do is
3092 * to initialize the `BaseSolver` and
override the two
virtual functions
3093 * (`make_mesh` and `system_matrix_local`).
3096 *
class Solver : public BaseSolver
3099 *
Solver() = delete;
3100 *
Solver(
const unsigned int p,
3101 *
const unsigned int mapping_degree,
3105 *
const std::string &file_name =
"data",
3109 *
virtual void make_mesh() override final;
3111 *
system_matrix_local(const CellIteratorPair &IP,
3112 *
AssemblyScratchData &scratch_data,
3113 *
AssemblyCopyData ©_data) override final;
3120 * Following is the implementation of the constructor.
3121 * We use the
second constructor of the `BaseSolver` class as
3122 * the solver is used at the
3123 * [fourth stage](@ref Step98_FourStages). By looking at the
3124 * [expressions](@ref Step98_Numerical_Recipe_B)
3125 *
for @f$A_{ij}@f$ and @f$b_i@f$ we can conclude that to compute them we need
3126 *
values of the shape
functions and the quadrature weights multiplied by the
3130 * numerically computed potential, @f$\vec{
A}@f$, on the right-hand side of the
3131 * [equation](@ref Step98_PDE_B). It is modeled by the
FE_Nedelec finite
3132 * elements. To compute the right-hand side, we need gradients of the shape
3137 *
Solver::Solver(
const unsigned int p,
3138 *
const unsigned int mapping_degree,
3142 *
const std::string &file_name,
3144 *
: BaseSolver(triangulation_ext,
3157 * At the fourth stage we
do not
load the mesh. We reuse the mesh loaded
3158 * at the
first stage. Consequently, we just need to distribute the dofs.
3161 *
void Solver::make_mesh()
3163 *
dof_handler.reinit(triangulation_rhs);
3164 *
dof_handler.distribute_dofs(fe);
3169 * The following function assembles a fraction of
3170 * [the system
matrix and the system right-hand side](@ref Step98_Numerical_Recipe_B)
3171 * related to a single cell. These fractions are
3172 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied to
3173 * `system_matrix` and `system_rhs` by
WorkStream.
3176 *
void Solver::system_matrix_local(
const CellIteratorPair &IP,
3177 *
AssemblyScratchData &scratch_data,
3178 *
AssemblyCopyData ©_data)
3182 *
copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
3184 *
scratch_data.dofs_per_cell);
3186 *
copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
3188 *
copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
3194 *
scratch_data.fe_values.reinit(cell);
3195 *
scratch_data.fe_values_rhs.reinit(cell_rhs);
3197 *
scratch_data.fe_values_rhs.get_function_gradients(
3198 *
scratch_data.dofs_rhs, scratch_data.vectors_vectors_list_rhs);
3200 *
for (
unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index)
3202 *
for (
unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i)
3204 *
for (
unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j)
3206 *
copy_data.cell_matrix(i, j) +=
3207 *
scratch_data.fe_values[se].value(i, q_index) *
3208 *
scratch_data.fe_values[se].value(j, q_index) *
3209 *
scratch_data.fe_values.JxW(q_index);
3212 *
copy_data.cell_rhs(i) +=
3213 *
(scratch_data.vectors_vectors_list_rhs[q_index][1][0] -
3214 *
scratch_data.vectors_vectors_list_rhs[q_index][0][1]) *
3215 *
scratch_data.fe_values[se].value(i, q_index) *
3216 *
scratch_data.fe_values.JxW(
3221 *
cell->get_dof_indices(copy_data.local_dof_indices);
3228 * <a name=
"step_98-Themainloop"></a>
3234 * function. The implementation of the
loop is straightforward -
3235 * we create and
run the four solvers
one-by-
one and
copy the relevant
3236 *
data into convergence tables.
3239 *
class MagneticProblem
3244 *
if (Settings::n_threads_max != 0)
3247 *
MainOutputTable table_T(2);
3248 *
MainOutputTable table_Jf(2);
3249 *
MainOutputTable table_A(2);
3250 *
MainOutputTable table_B(2);
3252 *
std::cout <<
"Solving for (p' = " << Settings::fe_degree + 1
3253 *
<<
"; p = " << Settings::fe_degree <<
"): " << std::flush;
3255 *
for (
unsigned int r = 1; r < 5; r++)
3257 *
table_T.add_value(
"r", r);
3258 *
table_T.add_value(
"p", Settings::fe_degree + 1);
3260 *
table_Jf.add_value(
"r", r);
3261 *
table_Jf.add_value(
"p", Settings::fe_degree);
3263 *
table_A.add_value(
"r", r);
3264 *
table_A.add_value(
"p", Settings::fe_degree);
3266 *
table_B.add_value(
"r", r);
3267 *
table_B.add_value(
"p", Settings::fe_degree);
3271 * Stage 1. Computing @f$T@f$.
3277 *
std::cout <<
"T " << std::flush;
3279 *
ExactSolutions::CurrentVectorPotential T_exact;
3281 *
SolverT::Solver
T(Settings::fe_degree + 1,
3283 *
Settings::mapping_degree,
3284 *
"T_p" + std::to_string(Settings::fe_degree + 1) +
3285 *
"_r" + std::to_string(r),
3290 *
table_T.add_value(
"ndofs",
T.get_n_dofs());
3291 *
table_T.add_value(
"ncells",
T.get_n_cells());
3292 *
table_T.add_value(
"L2",
T.get_L2_norm());
3296 * Stage 2. Computing @f$\vec{
J}_f@f$.
3302 *
std::cout <<
"Jf " << std::flush;
3304 *
ExactSolutions::FreeCurrentDensity Jf_exact;
3306 *
ProjectorHgradToHdiv::Solver Jf(Settings::fe_degree,
3307 *
Settings::mapping_degree,
3309 *
T.get_dof_handler(),
3312 *
std::to_string(Settings::fe_degree) +
3313 *
"_r" + std::to_string(r),
3318 *
table_Jf.add_value(
"ndofs", Jf.get_n_dofs());
3319 *
table_Jf.add_value(
"ncells", Jf.get_n_cells());
3320 *
table_Jf.add_value(
"L2", Jf.get_L2_norm());
3324 * Stage 3. Computing @f$\vec{
A}@f$.
3330 *
std::cout <<
"A " << std::flush;
3332 *
ExactSolutions::MagneticVectorPotential A_exact;
3334 *
SolverA::Solver
A(Settings::fe_degree,
3335 *
Settings::mapping_degree,
3337 *
T.get_dof_handler(),
3339 *
"A_p" + std::to_string(Settings::fe_degree) +
"_r" +
3340 *
std::to_string(r),
3345 *
table_A.add_value(
"ndofs",
A.get_n_dofs());
3346 *
table_A.add_value(
"ncells",
A.get_n_cells());
3347 *
table_A.add_value(
"L2",
A.get_L2_norm());
3351 * Stage 4. Computing @f$B@f$.
3357 *
std::cout <<
"B " << std::flush;
3359 *
ExactSolutions::MagneticField B_exact;
3361 *
ProjectorHcurlToL2::Solver B(Settings::fe_degree,
3362 *
Settings::mapping_degree,
3364 *
A.get_dof_handler(),
3367 *
std::to_string(Settings::fe_degree) +
3368 *
"_r" + std::to_string(r),
3372 *
table_B.add_value(
"ndofs", B.get_n_dofs());
3373 *
table_B.add_value(
"ncells", B.get_n_cells());
3374 *
table_B.add_value(
"L2", B.get_L2_norm());
3382 *
table_T.save(
"table_T_p" + std::to_string(Settings::fe_degree + 1));
3383 *
table_Jf.save(
"table_Jf_p" + std::to_string(Settings::fe_degree));
3384 *
table_A.save(
"table_A_p" + std::to_string(Settings::fe_degree));
3385 *
table_B.save(
"table_B_p" + std::to_string(Settings::fe_degree));
3386 *
std::cout << std::endl;
3390 *
int main(
int argc,
char **argv)
3396 * gmsh has global state, so
set that up in the normal way:
3402 *
MagneticProblem problem;
3405 *
catch (std::exception &exc)
3407 *
std::cerr << std::endl
3409 *
<<
"----------------------------------------------------"
3411 *
std::cerr <<
"Exception on processing: " << std::endl
3412 *
<< exc.what() << std::endl
3413 *
<<
"Aborting!" << std::endl
3414 *
<<
"----------------------------------------------------"
3421 *
std::cerr << std::endl
3423 *
<<
"----------------------------------------------------"
3425 *
std::cerr <<
"Unknown exception!" << std::endl
3426 *
<<
"Aborting!" << std::endl
3427 *
<<
"----------------------------------------------------"
3435<a name=
"step_98-Results"></a><h1>Results</h1>
3438The program generates the following output in the command line
interface by
3442Solving
for (p' = 1; p = 0):
T Jf
A B
T Jf
A B
T Jf
A B
T Jf
A B
3445The program assumes the finite elements of the lowermost degree, i.
e., @f$p
' = 1@f$
3446for the FE_Q finite elements and @f$p=0@f$ for other finite elements. To change the
3447degree of the finite elements, say @f$p' = 3@f$ and @f$p = 2@f$,
one needs to change
3448the setting `
Settings::fe_degree = 2` and rebuild the program. The degree of the
3449FE_Q finite elements will be computed automatically as @f$p
'= p + 1@f$.
3451The program dumps a number of files in the current directory. In the default
3452configuration these files are:
3453- `.vtu` files. They contain the computed vector fields. Recall that the spherical
3454 manifold and transfinite interpolation manifold are attached to many cell faces.
3455 Consequently, these cell faces are curved. Furthermore, the shape functions are
3456 mapped from the reference cell to the real mesh cells by the second-order mapping
3457 to accommodate the cells with curved faces. For these reasons, one needs to use a
3458 visualization software that can deal with curved faces and the higher-order mapping.
3459 A fresh version of ParaView is recommended. Visit did not have this feature
3460 at the time this tutorial was written (in early 2026). The
3461 <a href="https://github.com/dealii/dealii/wiki/Notes-on-visualizing-high-order-output"> Notes on visualizing high order output</a>
3462 provide more information on this topic.
3463- `.tex` files. These files contain the convergence tables.
3465The following provides examples of the convergence tables simulated with the
3466default settings for three different degrees of the finite elements,
3467@f$p = 0, 1, 2@f$ (recall that @f$p' = p + 1@f$).
3470<caption>Convergence table @f$T@f$.</caption>
3476 <th>@f$\|e\|_{L^2}@f$</th>
3477 <th>@f$\alpha_{L^2}@f$</th>
3578<caption>Convergence table @f$\vec{J}_f@f$.</caption>
3584 <th>@f$\|e\|_{L^2}@f$</th>
3585 <th>@f$\alpha_{L^2}@f$</th>
3686<caption> Convergence table @f$B@f$.</caption>
3692 <th>@f$\|e\|_{L^2}@f$</th>
3693 <th>@f$\alpha_{L^2}@f$</th>
3793The following notations were used in the headers of the tables:
3795- p', p - the degree of the finite elements.
3797- r - the mesh refinement parameter, i.
e., the number of global mesh
3800- cells - the total number of active cells.
3802- dofs - the number of degrees of freedom.
3804-@f$\|
e\|_{
L^2}@f$ - the @f$L^2@f$ error
norm.
3806-@f$\alpha_{
L^2}@f$ - the order of convergence of the @f$L^2@f$ error
norm.
3808Let us contemplate these convergence tables
for a brief moment. The
first table
3809illustrates convergence of the numerically computed current vector potential,
3810@f$T@f$. In this particular case we can expect the order of the convergence rate of
3811@f$\alpha_{
L^2} \le p
' + 1@f$ (See also <a href="https://www.math.colostate.edu/~bangerth/videos.676.3.95.html">video lecture 3.95</a>.) If the finite
3812elements of the lowermost degree are used, @f$p'=1@f$, the order of
3813convergence rate is at the upper boundary of the expected
values,
3814@f$\alpha_{
L^2} \approx 2.0@f$. If the finite elements of the
second degree are
3815used, @f$p
'=2@f$, the order of the convergence rate is higher than expected,
3816@f$\alpha_{L^2} \approx 4.0@f$. This means that the @f$L^2@f$ error norm converges to
3817zero at a rate higher than theoretically possible. Most likely this is due to
3818the fact that the current vector potential has a particularly simple form,
3819[Equation 3](@ref Step98_Equation_3_T). It is either constant or changes as
3820the second-order monomial,
3824The second order polynomial approximates this behavior exactly. That is to
3825say, in this particular case we have a lucky situation in which the shape
3826functions can approximate the numerical solution exactly within each mesh
3827cell. This, most likely, explains the extra rapid convergence of the @f$L^2@f$
3828norm. Note also that the error norms, @f$\|e\|_{L^2}@f$, at @f$p'=2@f$ and @f$p
'=3@f$
3829are the same for the same values of the mesh refinement parameter, @f$r@f$. This
3830is, most likely, due to the fact that the second-order shape functions,
3831@f$p'=2@f$,
approximate the solution exactly within each mesh cell and the
3832third-order monomials of the shape function at @f$p
'=3@f$ have absolutely
3833nothing to contribute to the quality of approximation.
3835The second table illustrates convergence of the numerically computed
3836free-current density, @f$\vec{J}_f@f$. The free-current density has been
3837computed as a derivative of current vector potential,
3839\vec{J}_f = \vec{\nabla}\overset{V} {\times} T.
3841The derivative reduces the order of the convergence rate by one. Most
3842likely in this particular case the error made in computing @f$\vec{J}_f@f$
3843is defined by the error made in computing @f$T@f$. For this reason, the order
3844of convergence in the second table equals the order of convergence in
3845the first table minus one.
3847Due to the implicit gauge we cannot observe the convergence of the magnetic
3848vector potential, @f$\vec{A}@f$. Instead, we can observe the convergence on the
3849magnetic field, @f$B@f$, given in the third table. It follows the same pattern:
3850The rate of convergence at the lower degrees of the finite elements is at
3851the best expected value, @f$\alpha_{L^2} = p + 1@f$; the rate of convergence at
3852the higher degrees of the finite elements is better than expected. The extra
3853rapid convergence at the higher degrees of the finite elements is explained
3854by the relatively simple form of the field being approximated.
3856The figures below illustrate the current vector potential, @f$T@f$, the free-current
3857density, @f$\vec{J}_f@f$, and magnetic field, @f$B@f$, computed with the following
3858settings: @f$p = 2@f$ and @f$r = 4@f$. Visual inspection of the magnetic vector
3859potential, @f$\vec{A}@f$, is not very informative as its conservative portion is
3864 <img src="https://dealii.org/images/steps/developer/step-98-result-T.svg"
3865 alt="The result - current vector potential" height="531">
3871 <img src="https://dealii.org/images/steps/developer/step-98-result-Jf.svg"
3872 alt="The result - free current density" height="531">
3878 <img src="https://dealii.org/images/steps/developer/step-98-result-B.svg"
3879 alt="The result - magnetic field" height="531">
3883The images above suggest that the computed fields do not exhibit any irregular
3884behavior (the first image in the table below illustrates how irregular
3885behavior can look like). The fields on these images closely resemble the
3886corresponding closed-form analytical expressions given in the introduction.
3887From the first glans the convergence tables above may appear somewhat strange
3888due to extra rapid convergence at the higher degrees of the finite elements.
3889One, however, can argue that the extra rapid convergence can be explained by
3890the simple polynomial form of the fields being approximated. One thing is
3891certain - the convergence rates presented in these tables are at the best
3892theoretically expected values or better.
3894@anchor Step98_PossibilitiesForExtensions
3895<a name="step_98-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
3897Let us consider the two-dimensional curl-curl partial differential equation
3900\vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg)
3901+ \eta^2 \vec{A} = \vec{\nabla}\overset{V}{\times} T.
3903We can compute the divergence of this expression and rearrange the terms as the
3906\vec{\nabla}\cdot \vec{A} = \frac{1}{\eta^2} \vec{\nabla}\cdot
3907\bigg[\vec{\nabla}\overset{V}{\times} T -
3908\vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu}
3909\vec{\nabla}\overset{S}{\times}\vec{A}\bigg)\bigg].
3911The divergence of the vector curl equals zero, see the introduction. Therefore,
3912the right-hand side of the last equation equals zero,
3914\vec{\nabla}\cdot \vec{A} = 0.
3916That is to say, the @f$\eta^2@f$ gauging term can be considered to be the Coulomb
3917gauge at least in theory. In practice the situation is a bit more complicated.
3918What the @f$\eta^2@f$ gauging term does depends on the value of @f$\eta^2@f$. The table
3919below attempts to express this very point.
3922 <img src="https://dealii.org/images/steps/developer/step-98-eta.svg"
3923 alt="The result - eta-squared table" height="1062">
3927This table presents four simulations for four different values of @f$\eta^2@f$. In
3928all four simulations the degree of the finite elements and the mesh refinement
3929parameter were @f$p=2@f$ and @f$r=1@f$, respectively. The setting
3930`Settings::project_exact_solution` was set to `true`.
3932At the setting @f$\eta^2=\dfrac{10^{-10}}{\mu_0}@f$ the computed magnetic vector
3933potential (the orange curve on the @f$|\vec{A}(x,0)|@f$ plot) looks exactly the same
3934as the projected exact solution gauged by the Coulomb gauge (the blue curve). At
3935this value of @f$\eta^2@f$ the term @f$\eta^2\vec{A}@f$ acts like the Coulomb gauge.
3936At the default setting, @f$\eta^2=0@f$, the solution, @f$\vec{A}@f$, is contaminated by
3937an unknown conservative vector field. The computed magnetic field, @f$B@f$, at this
3938setting looks exactly as the corresponding exact solution because the
3939conservative portion of @f$\vec{A}@f$ is filtered out by the process of computing
3940the magnetic field, @f$B = \vec{\nabla}\overset{S}{\times}\vec{A}@f$. At the setting
3941@f$\eta^2=\dfrac{10^{-6}}{\mu_0}@f$ the solenoidal part of the computed magnetic
3942vector potential deviates from the exact solution. We can deduce this from the
3943fact that the computed magnetic
3944field deviates from the exact expression for the magnetic field on the @f$B(x,0)@f$
3945plot. The error in the solenoidal part of the computed @f$\vec{A}@f$ is due to
3946the fact that introduction of the gauging term, @f$\eta^2\vec{A}@f$, modifies the
3948curl-curl equation. So, strictly speaking, we are solving a different partial
3949differential equation. If @f$\eta^2@f$ is small, the difference between the initial
3950curl-curl equation and the curl-curl equation modified by adding the gauging
3951term is negligible. Consequently, the error in the solenoidal part of the
3952computed @f$\vec{A}@f$ is negligible as well. Evidently, in this particular case
3953"small" means @f$\eta^2 \ll \dfrac{10^{-6}}{\mu_0}@f$.
3955Normally, we are interested in measurable fields such as magnetic field and
3956consider the magnetic vector potential as a useful tool for computing measurable
3957fields. We can tolerate the presence of an unknown conservative component in the
3958magnetic vector potential, i.e., the situation illustrated by the first two rows
3959in the table above. In such disposition we need to keep @f$\eta^2@f$ as small as
3960possible, i.e., as far away as possible from the situation shown in last row of
3961the table. In this tutorial program we set @f$\eta^2@f$ to zero and increase it just
3962a bit in if the conjugate gradient algorithm cannot converge.
3964Suppose for a moment that we would like to have the solution to the curl-curl
3965equation in terms of a purely solenoidal magnetic vector potential, @f$\vec{A}@f$,
3966that is, a solution gauged by the Coulomb gauge. To get such a solution we need
3967to tweak the @f$\eta^2@f$ parameter. By contemplating
3968the table above one can hypothesize that there exists an optimal value of the
3969gauging parameter, @f$\eta^2_\text{opt}@f$, at which the @f$L^2@f$ error norm computed for
3970@f$\vec{A}@f$ is minimal. The optimal value should be somewhere in between
3971@f$\eta^2=\dfrac{10^{-12}}{\mu_0}@f$ and @f$\eta^2=\dfrac{10^{-6}}{\mu_0}@f$, based on the
3972experiments above. Try to
3973verify this hypothesis by finding the exact value of @f$\eta^2_\text{opt}@f$.
3975In this instance we have a close-form analytical expression of the exact
3976solution, i.e., the expression of @f$\vec{A}@f$ given in the introduction. In a
3977real-life simulation there is no expression of the exact solution. Try to
3978think of a method of blind (meaning without the exact solution) estimation of
3979@f$\eta^2_\text{opt}@f$. Try to implement and test your ideas.
3981Adding the gauging term, @f$\eta^2\vec{A}@f$, to the curl-curl equation converts
3982a positive semidefinite system matrix into a positive definite matrix. Strictly
3983speaking, adding the gauging term modifies the initial curl-curl equation. For
3984this reason, it is important to keep @f$\eta^2@f$ small so the solution is not
3985afflicted by the error induced by adding the gauging term. The definition of
3986"small" here is a bit fuzzy. In absence of the exact solution setting @f$\eta^2@f$
3987at the acceptable level of the error in the solenoidal component of @f$\vec{A}@f$
3988is difficult. It is better to discard the gauging term and use a more
3989sophisticated linear solver. The hypre AMS @cite hypre1998b can solve the
3990systems of linear equations yielded by the curl-curl equation without the
3991@f$\eta^2@f$ gauging term, (@f$\beta=0@f$ in @cite hypre1998b). Try to implement the
3995<a name="step_98-PlainProg"></a>
3996<h1> The plain program</h1>
3997@include "step-98.cc"
* * for(const auto &cell :triangulation.active_cell_iterators())
* * int main(int argc, char **argv)
* x_component_mask set(0, true)
* * * struct InterferenceTaperTransform *
ConvergenceTable()=default
void evaluate_convergence_rates(const std::string &data_column_key, const std::string &reference_column_key, const RateMode rate_mode, const unsigned int dim=2)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< RangeNumberType > > &values) const
void attach_triangulation(Triangulation< dim, spacedim > &tria)
static void set_thread_limit(const unsigned int max_threads=numbers::invalid_unsigned_int)
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
void write_tex(std::ostream &file, const bool with_header=true) const
void set_column_order(const std::vector< std::string > &new_order)
void set_tex_caption(const std::string &key, const std::string &tex_caption)
void set_scientific(const std::string &key, const bool scientific)
void set_precision(const std::string &key, const unsigned int precision)
void initialize(const Triangulation< dim, spacedim > &triangulation)
#define Assert(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
@ update_default
No update.
std::vector< index_type > data
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
DataComponentInterpretation
@ component_is_part_of_vector
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
* * * ScaleZFunction< dim, Number, components >::ScaleZFunction * component(component)
* * if(update_pressure &update_flags) * compute_pressure(constitutive_request
* * * * * * TimeRateUpdateFlags TimeRateRequest< ValueType, dim, Number > get_update_flags() const
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, ShapeDataMemorySpace > shape_data, const ViewTypeIn in, ViewTypeOut out)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
void load(Archive &ar, ::std_cxx26::inplace_vector< T, N > &vec, const unsigned int)
void save(Archive &ar, const ::std_cxx26::inplace_vector< T, N > &vec, const unsigned int)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
unsigned int get_degree(const std::vector< typename BarycentricPolynomials< dim >::PolyType > &polys)
constexpr types::manifold_id flat_manifold_id
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
bool write_higher_order_cells
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
1.9.8