1517 *
const unsigned int mapping_degree = 2;
1518 *
const unsigned int fe_degree = 0;
1520 *
const double eta_squared = 0.0;
1522 *
const unsigned int n_threads_max = 0;
1524 *
const bool project_exact_solution =
false;
1530 * <a name=
"step_98-Convergencetable"></a>
1531 * <h3>Convergence table</h3>
1532 * The following
class describes a convergence table. The convergence tables are
1533 * saved on disk in TeX format.
1539 *
MainOutputTable() =
delete;
1541 *
MainOutputTable(
const unsigned int dim)
1546 *
void save(
const std::string &file_name)
1565 *
std::ofstream ofs(file_name +
".tex");
1570 *
const unsigned int dim;
1576 * <a name=
"step_98-Equations"></a>
1577 * <h3>Equations</h3>
1578 * The following
namespace contains closed-form analytical expressions
1579 *
for @f$T@f$, @f$\vec{
J}_f@f$, @f$\vec{
A}@f$, @f$B@f$, mentioned in the introduction to
this
1583 *
namespace ExactSolutions
1588 * The following function describes the free-current density, @f$\vec{
J}_f@f$,
1589 *
inside the current region. The current density in
this tutorial is
1590 * implemented by `ExactSolutions::FreeCurrentDensity`
class and by
1591 * `SolverT::Solver::free_current_density` member function . There is a
1592 * subtle difference in how these two implementation compute the free-current
1593 * density. Both classes, however, utilize the same expression
for the
1594 * free-current density. This function describes the expression.
1605 * The following
class implements the closed-form analytical
1606 * [expression](@ref Step98_Equation_1_Jf)
for the free-current density,
1607 * @f$\vec{
J}_f@f$, in the entire domain. The free-current density is computed
1608 * purely on the basis of the spatial coordinates of the field
point. In go
1609 * coordinates, out comes the free-current density. The information on the
1610 * material ID of the mesh cells and any other information on the mesh is
1611 * ignored. This function is used
for computing @f$L_2@f$ error norms and
for
1612 * computing the projected exact solution. The @f$\vec{
J}_f@f$ on the right-hand
1613 * side of the div-grad equation is implemented by the member function
1614 * `SolverT::Solver::free_current_density`
1617 *
class FreeCurrentDensity :
public Function<2>
1620 *
FreeCurrentDensity()
1629 *
ExcDimensionMismatch(
values.size(), p.size()));
1631 *
for (
unsigned int i = 0; i <
values.size(); i++)
1633 *
const double r = p[i].norm();
1635 *
if ((r >= Settings::a2) && (r <= Settings::b2))
1638 *
volume_free_current_density(p[i], Settings::K0);
1651 * The following
class implements the closed-form analytical
1652 * [expression](@ref Step98_Equation_4_A)
1653 *
for the magnetic vector potential, @f$\vec{
A}@f$.
1656 *
class MagneticVectorPotential :
public Function<2>
1659 *
MagneticVectorPotential()
1668 *
ExcDimensionMismatch(
values.size(), p.size()));
1670 *
for (
unsigned int i = 0; i <
values.size(); i++)
1672 *
const double r = p[i].norm();
1679 *
A = (mu_1 * K0 / 4.0) * (b2 * b2 - a2 * a2);
1683 *
A = (mu_0 * K0 / 4.0) * (b2 * b2 - a2 * a2) *
1684 *
(r + b1 * b1 * (mu_r - 1.0) / r) / r;
1688 *
A = (mu_0 * K0 / 2.0) *
1689 *
(b2 * b2 * r / 2.0 -
std::pow(r, 3) / 4.0 +
1690 *
(a2 * (b2 * b2 - a2 * a2) *
1691 *
(a2 + b1 * b1 * (mu_r - 1.0) / a2) / 2.0 -
1698 *
A = (mu_0 * K0 / 2.0) * b2 *
1700 *
(a2 * (b2 * b2 - a2 * a2) *
1701 *
(a2 + b1 * b1 * (mu_r - 1.0) / a2) / 2.0 -
1707 *
values[i][0] = -
A * p[i][1];
1715 * The following
class implements the closed-form analytical
1716 * [expression](@ref Step98_Equation_3_T)
for
1717 * the current vector potential, @f$T@f$.
1720 *
class CurrentVectorPotential : public
Function<2>
1723 *
CurrentVectorPotential()
1736 *
ExcMessage(
"This line is to avoid compiler warnings."));
1738 *
for (
unsigned int i = 0; i <
values.size(); i++)
1740 *
const double r = p[i].norm();
1745 *
values[i] = K0 * (b2 * b2 - a2 * a2) / 2.0;
1749 *
values[i] = -K0 * (r * r - b2 * b2) / 2.0;
1761 * The following
class implements the closed-form analytical
1762 * [expression](@ref Step98_Equation_2_B)
1763 *
for the magnetic field, @f$B@f$.
1766 *
class MagneticField : public
Function<2>
1782 *
ExcMessage(
"This line is to avoid compiler warnings."));
1784 *
for (
unsigned int i = 0; i <
values.size(); i++)
1786 *
const double r = p[i].norm();
1791 *
values[i] = mu_1 * K0 * (b2 * b2 - a2 * a2) / 2.0;
1795 *
values[i] = mu_0 * K0 * (b2 * b2 - a2 * a2) / 2.0;
1799 *
values[i] = mu_0 * K0 * (b2 * b2 - r * r) / 2.0;
1814 * <a name=
"step_98-BaseSolver"></a>
1815 * <h3>Base Solver</h3>
1819 * As discussed above, the following
namespace aggregates the code common to all
1820 * four solvers used in the tutorial. All four solvers are derived from the
1821 * `BaseSolver`
class.
1824 *
namespace BaseClasses
1828 * The computation of the integrals of the functionals is delegated to the
1829 * derived classes. The function that computes the integrals,
1830 * `BaseSolver::system_matrix_local`, is
virtual and must be overridden by
1831 * the derived classes. However, the objects of the
types FEValues are
1832 * initialized at the
level of the `BaseSolver`
class together with
1833 * `AssemblyScratchData`. For the initialization to work properly, the
1834 * `BaseSolver`
class must know which cell
data to compute
for a particular
1835 * implementation of the solver down the hierarchy. This information is
1836 * communicated to the `BaseSolver` by passing an argument of the type
1837 * `UpdateFlagsCollection` to the constructor. In all solvers, with exception
1838 * of `SolverT`, we have two
types of finite elements. One type of finite
1839 * elements models the solution to the partial differential equation. Another
1840 * models the physical quantity on the right-hand side of the partial
1841 * differential equation. The `solution_update_flags`
data member below
1842 * contains the flags
for updating the
values of the finite elements that
1843 * model the solution. The `rhs_update_flags`
data member contains the flags
1844 *
for updating the
values of the finite elements that model the physical
1845 * quantity on the right-hand side of the partial differential equation.
1848 *
struct UpdateFlagsCollection
1856 * Each iteration of the program consists of
1857 * [four stages](@ref Step98_FourStages). Each stage utilizes one solver. All
1858 * four solvers used in
this tutorial are derived from the following
class.
1859 * The solver used in the
first stage loads the mesh and refines it
if
1860 * necessary. The solvers in the other three stages reuse the mesh prepared at
1861 * the
first stage. Furthermore, the solver in the
first stage expects a
1862 * closed-form analytical expression on the right-hand side of the partial
1863 * differential equation. Each solver in the other three stages expects a
1864 * potential computed at one of the preceding stages, i.e., a field in a form
1865 * of linear superposition of the shape
functions. At the top of the following
1866 *
class, we declare two constructors. The
first constructor must be used
for
1867 * constructing solver
for the
first stage. The
second constructor must be
1868 * used
for constructing the solvers
for the
second, third, and fourth stages.
1869 * The
second constructor has three extra arguments, `triangulation_rhs`,
1870 * `dof_handler_rhs`, and `solution_rhs`
for accommodating the mesh and the
1871 * potential computed at one of the preceding stages.
1875 * Following the constructors are the eight functions that implement various
1876 * steps typical
for every solver. The function `run` aggregates these steps.
1877 * This arrangement of functions is quite standard in deal.II, see @ref step_3
"step-3",
1878 *
for instance. Normally, the `setup` function begins by distributing the
1879 * dofs. The problem is: the type of the finite elements is not known in the
1880 * `BaseSolver` class. It is specified in the derived classes. Therefore, it
1881 * could be reasonable to make the setup function virtual as well. Instead,
1882 * we move the dof distribution code to the
end of the `make_mesh` function
1883 * which is, in fact, virtual and must be overridden in the derived classes
1888 * After that, we declare six get
functions that provide access to protected
1890 * `MagneticProblem::
run` function. The
data provided by these get-
functions
1891 * is used to fill the convergence tables and to pass the references to the
1892 * triangulation, the dof handler, and the dofs between solvers.
1898 *
BaseSolver(
const unsigned
int mapping_degree,
1899 *
const UpdateFlagsCollection update_flags_collection,
1906 *
const unsigned int stage,
1907 *
const unsigned int mapping_degree,
1908 *
const UpdateFlagsCollection update_flags_collection,
1909 *
const std::string &file_name,
1912 *
virtual void make_mesh() = 0;
1919 *
void output_results()
const;
1920 *
void compute_error_norms();
1921 *
void project_exact_solution_fcn();
1926 *
double get_L2_norm()
const;
1927 *
unsigned int get_n_cells()
const;
1935 * We
begin the `
protected` section of the `BaseSolver`
class by declaring
1936 * three
data members that store the input from
one of the preceding
1937 * solvers. If there is no preceding solver and these
data is not provided,
1938 * i.e.,
first of the two constructors above has been used, these three
data
1939 * members
point to `triangulation`, `dof_handler`, and `solution`, of the
1944 * Following are the three
data members that store the triangulation, dof
1945 * handler, and dofs vector of the current solver. In the
case of the
1946 *
first-stage solver, the `triangulation`
data member stores the loaded and
1947 * refined mesh. This
data member is not used in the
case of the solvers of
1948 * the
second, third, and fourth stages. The `dof_handler` and `solution`
1949 * represent the result of the solver, i.e., the computed potential or
1954 * Next, we declare four
data members that describe the system of linear
1955 * equations to be solved by the linear solver. The components of the system
1956 *
matrix and that of the right-hand side are computed by the `
assemble`
1957 * function. The
affine constraints are used to
apply the Dirichlet boundary
1958 * conditions and to distribute the local (cell specific) system
matrix and
1959 * right-hand side to `system_matrix` and `system_rhs`. The are no hanging
1960 * nodes and hanging node constraints in
this program. The last
data member
1961 * in
this block describes the
dynamic sparsity pattern. The tutorial @ref step_2
"step-2"
1962 * discusses the rationale behind the
dynamic sparsity pattern.
1966 * The next block contains two
data members related to the exact solution.
1967 * The
data member `exact_solution` points to the closed-form analytical
1968 * solution the solver attempts to compute.
1969 * If `Settings::project_exact_solution=
true`, the exact solution is
1970 * projected onto a proper function space and the
data member
1971 * `projected_exact_solution` is populated by the dofs of the projected
1972 * exact solution. The corresponding dof handler is `dof_handler`. Together
1973 * `dof_handler` and `projected_exact_solution` constitute the field
1974 * function which describes the exact solution. It is saved into the `.vtu`
1975 * file next to the solution and the @f$L_2@f$ error
norm.
1979 * The next four
data members simply store the
data supplied as arguments
1980 * to the constructor. The `stage`
data member stores the number of the
1981 * current stage. The `mapping_degree`
data member contains the degree of
1982 * mapping from the
reference cell to a mesh cell and back. The
1983 * `update_flags_collection` contains the information on which finite
1984 * element
values must be computed
for each cell. The names of the output
1985 * files are derived by appending strings to `file_name`.
1989 * The next block contains two
data members computed by the function
1990 * `BaseSolver::compute_error_norms`. The `L2_per_cell`
data member
1991 * contains
one value of the @f$L^2@f$ error
norm per mesh cell. It is saved
1992 * into the `.vtu` file next to the solution. The `
L2_norm`
data member
1993 * contains
one value of the @f$L^2@f$ error
norm per mesh. It is reported in
1994 * the convergence table.
2014 *
const unsigned int stage;
2015 *
const unsigned int mapping_degree;
2016 *
const UpdateFlagsCollection update_flags_collection;
2017 *
const std::string file_name;
2024 * The program utilizes the
WorkStream technology. The @ref step_9
"step-9" tutorial
2025 * does a much better job of explaining the workings of
WorkStream.
2026 * Reading the @ref workstream_paper
"WorkStream paper" is recommended.
2027 * In very simple terms, the workings of the
WorkStream can be envisioned as
2028 * the following. Let us assume we have a task of computing components of
2029 * the system
matrix, @f$A_{ij}@f$, and the right-hand side vector, @f$b_i@f$.
2030 * Simply put, we need to fill in the
matrix `system_matrix` and vector
2031 * `system_rhs`. This is a big task as the number of dofs is large. The
2032 * idea is to
split the big task on a number of small tasks and feed them
2033 * to multiple threads to speed up the calculation process. Each small
2034 * task consists of computing the contributions of a single mesh cell to
2035 * `system_matrix` and `system_rhs`. These contributions are stored
2036 * temporary in `
cell_matrix` and `cell_rhs`
for each cell. These
2037 * contributions are then copied to `system_matrix` and `system_rhs`.
2038 *
WorkStream creates and schedules the small tasks and takes care
2039 * of copying `
cell_matrix` and `cell_rhs` to `system_matrix` and
2040 * `system_rhs`. The rest of the declarations in the `
protected` section
2041 * of the `BaseSolver`
class help to communicate to
WorkStream information
2042 * which is necessary
for its operation.
2046 * First, the `CellIteratorPair` type is declared in the block of code
2047 * below. The solvers in the
second, third, and fourth stages use two dof
2048 * handlers, `dof_handler` and `dof_handler_rhs`. The
WorkStream needs to
2049 * walk through the two dof handlers synchronously. For
this purpose we
2050 * pair two active cell iterators (one from `dof_handler`, another from
2051 * `dof_handler_rhs`). For that we need the `CellIteratorPair` type. The
2052 * solver at the
first stage (a solver constructed by invoking the
first
2053 * constructor above) uses only
one dof handler. In
this case the
2054 * constructor makes `dof_handler_rhs` to
reference `dof_handler`. In
2055 * effect, both iterators of the tuple will iterate the same dof handler.
2056 * In the
case of the
first-stage solver we will use only the
first
2061 * Next, we declare the `AssemblyScratchData` type. An
object of
this type
2062 * contains the relevant information on the current mesh cell which is used
2063 * as an input
for computing the components of the system
matrix and the
2064 * right-hand side.
WorkStream creates an
object of
this type and passes
2065 * it to the function `system_matrix_local` which, in turn, computes
2070 * Next, the type `AssemblyCopyData` is declared. An objects of
this type
2071 * contains the cell specific contributions to the system
matrix and the
2072 * right-hand side. The
WorkStream creates
object of
this type and passes
2073 * it to function `system_matrix_local` along with an
object of the type
2074 * `AssemblyScratchData`. The function `system_matrix_local`, in turn,
2075 * takes input
data form the
object of the type `AssemblyScratchData`,
2076 * computes the relevant integrals and places the result into the
object
2077 * of the type `AssemblyCopyData`. The
WorkStream copies the content of
2078 * the `AssemblyCopyData` object, `
cell_matrix` and `cell_rhs`, into
2079 * `system_matrix` and `system_rhs`. The
data member
2080 * `AssemblyCopyData::local_dof_indices` contains the indices of global
2081 * components to which cell-specific
data must be copied.
2085 *
WorkStream calls the function `system_matrix_local` to compute the
2086 * cell-specific components. Likewise,
WorkStream uses calls to
2087 * `copy_local_to_global` to
copy the cell-specific
data into the system
2089 * `copy_local_to_global` are declared last.
2092 *
using IteratorTuple =
2093 *
std::tuple<typename DoFHandler<2>::active_cell_iterator,
2098 *
struct AssemblyScratchData
2103 *
const unsigned int mapping_degree,
2104 *
const UpdateFlagsCollection update_flags_collection,
2105 *
const unsigned int stage);
2107 *
AssemblyScratchData(
const AssemblyScratchData &scratch_data);
2114 *
const unsigned int dofs_per_cell;
2115 *
const unsigned int n_q_points;
2117 *
std::vector<double> permeability_list;
2118 *
std::vector<double> values_list_rhs;
2119 *
std::vector<Tensor<1, 2>> vectors_list_rhs;
2120 *
std::vector<std::vector<Tensor<1, 2>>> vectors_vectors_list_rhs;
2126 *
struct AssemblyCopyData
2130 *
std::vector<types::global_dof_index> local_dof_indices;
2133 *
virtual void system_matrix_local(
const CellIteratorPair &IP,
2134 *
AssemblyScratchData &scratch_data,
2135 *
AssemblyCopyData ©_data) = 0;
2137 *
void copy_local_to_global(
const AssemblyCopyData ©_data);
2142 * The following are the implementations of the two constructors
2143 * of the `BaseSolver`
class.
2146 *
BaseSolver::BaseSolver(
const unsigned int mapping_degree,
2147 *
const UpdateFlagsCollection update_flags_collection,
2148 *
const std::string &file_name,
2150 *
: triangulation_rhs(triangulation)
2151 *
, dof_handler_rhs(dof_handler)
2152 *
, solution_rhs(solution)
2153 *
, exact_solution(exact_solution)
2155 *
, mapping_degree(mapping_degree)
2156 *
, update_flags_collection(update_flags_collection)
2157 *
, file_name(file_name)
2163 *
const unsigned int stage,
2164 *
const unsigned int mapping_degree,
2165 *
const UpdateFlagsCollection update_flags_collection,
2166 *
const std::string &file_name,
2168 *
: triangulation_rhs(triangulation_rhs)
2169 *
, dof_handler_rhs(dof_handler_rhs)
2170 *
, solution_rhs(solution_rhs)
2171 *
, exact_solution(exact_solution)
2173 *
, mapping_degree(mapping_degree)
2174 *
, update_flags_collection(update_flags_collection)
2175 *
, file_name(file_name)
2180 * The following function applies the Dirichlet boundary condition, sets
2181 * up a sparsity pattern, and initializes the vectors and matrices. It is
2182 * common
for the setup function to distribute the dofs. The type of the
2183 * finite elements, however, is not known at the
level of `BaseSolver`
class.
2184 * The type of the finite elements is chosen in the derived classes. For
2185 *
this reason, the task of distributing the dofs is shifted to the
end of
2186 * the `make_mesh` function which is a
virtual function.
2190 * The Dirichlet boundary condition must be enforced only in the
first
2191 * [stage](@ref Step98_FourStages)
2192 * where the div-grad equation is solved
for the current vector potential,
T.
2193 * For
this reason we have `
if (stage == 1)` filter in the beginning of the
2194 * function. The boundary
value problem
for the curl-curl equation utilizes
2195 * the Neumann boundary condition. It is a natural boundary condition. It is
2196 * enforced by minimization of the functional. The two projectors,
2197 * @f$T \rightarrow \vec{
J}_f@f$ and @f$\vec{
A} \rightarrow B@f$, use no boundary
2201 *
void BaseSolver::setup()
2203 *
constraints.clear();
2210 *
Settings::outer_boundary_id,
2214 *
constraints.close();
2218 *
sparsity_pattern.copy_from(dsp);
2220 *
system_matrix.reinit(sparsity_pattern);
2221 *
solution.reinit(dof_handler.n_dofs());
2222 *
system_rhs.reinit(dof_handler.n_dofs());
2224 *
if (Settings::project_exact_solution && exact_solution)
2225 *
projected_exact_solution.reinit(dof_handler.n_dofs());
2227 *
if (exact_solution)
2228 *
L2_per_cell.reinit(triangulation.n_active_cells());
2233 * Formally, the following function assembles the system of linear equations.
2234 * In reality, however, it just spells all the magic words to get the
2235 *
WorkStream going. The interesting part, i.e., computing the components of
2236 * the system
matrix and the right-hand side, happens in the function
2237 * `system_matrix_local` which is a
virtual function. That is to say, the
2238 * recipe
for the functional is implemented in the derived
class by overriding
2239 * the
virtual function `system_matrix_local`.
2242 *
void BaseSolver::assemble()
2245 *
dof_handler_rhs.begin_active()}),
2247 *
{dof_handler.end(), dof_handler_rhs.end()}),
2249 *
&BaseSolver::system_matrix_local,
2250 *
&BaseSolver::copy_local_to_global,
2251 *
AssemblyScratchData(dof_handler,
2255 *
update_flags_collection,
2257 *
AssemblyCopyData());
2262 * The following are the implementation of the constructors of the
2263 * `AssemblyScratchData`. The
first constructor initializes the scratch
2264 *
data from the input parameters. The
second - from another
object of
2265 * the same type, i.e., a
copy constructor.
2268 *
BaseSolver::AssemblyScratchData::AssemblyScratchData(
2272 *
const unsigned int mapping_degree,
2273 *
const UpdateFlagsCollection update_flags_collection,
2274 *
const unsigned int stage)
2275 *
: mapping(mapping_degree)
2276 *
, fe_values(mapping,
2277 *
dof_handler.get_fe(),
2278 *
QGauss<2>((stage == 1) ? (dof_handler.get_fe().degree + 1) :
2279 *
(dof_handler.get_fe().degree + 2)),
2280 *
update_flags_collection.solution_update_flags)
2281 *
, fe_values_rhs(mapping,
2282 *
dof_handler_rhs.get_fe(),
2283 *
QGauss<2>((stage == 1) ? (dof_handler.get_fe().degree + 1) :
2284 *
(dof_handler.get_fe().degree + 2)),
2285 *
update_flags_collection.rhs_update_flags)
2286 *
, dofs_per_cell(fe_values.dofs_per_cell)
2287 *
, n_q_points(fe_values.get_quadrature().
size())
2288 *
, permeability_list(n_q_points)
2289 *
, values_list_rhs(n_q_points)
2290 *
, vectors_list_rhs(n_q_points)
2291 *
, vectors_vectors_list_rhs(n_q_points,
std::vector<
Tensor<1, 2>>(2))
2292 *
, dof_handler_rhs(dof_handler_rhs)
2293 *
, dofs_rhs(dofs_rhs)
2296 *
BaseSolver::AssemblyScratchData::AssemblyScratchData(
2297 *
const AssemblyScratchData &scratch_data)
2298 *
: mapping(scratch_data.mapping.
get_degree())
2299 *
, fe_values(mapping,
2300 *
scratch_data.fe_values.get_fe(),
2301 *
scratch_data.fe_values.get_quadrature(),
2303 *
, fe_values_rhs(mapping,
2304 *
scratch_data.fe_values_rhs.get_fe(),
2305 *
scratch_data.fe_values_rhs.get_quadrature(),
2307 *
, dofs_per_cell(fe_values.dofs_per_cell)
2308 *
, n_q_points(fe_values.get_quadrature().
size())
2309 *
, permeability_list(n_q_points)
2310 *
, values_list_rhs(n_q_points)
2311 *
, vectors_list_rhs(n_q_points)
2312 *
, vectors_vectors_list_rhs(n_q_points,
std::vector<
Tensor<1, 2>>(2))
2313 *
, dof_handler_rhs(scratch_data.dof_handler_rhs)
2314 *
, dofs_rhs(scratch_data.dofs_rhs)
2319 * The following function copies the components of a cell
matrix and a cell
2320 * right-hand side into the system
matrix, @f$A_{ij}@f$, and the system right-hand
2324 *
void BaseSolver::copy_local_to_global(
const AssemblyCopyData ©_data)
2326 *
constraints.distribute_local_to_global(copy_data.
cell_matrix,
2327 *
copy_data.cell_rhs,
2328 *
copy_data.local_dof_indices,
2335 * The following function solves the system of linear equations.
2336 * The stopping condition
for the iteration algorithm is
2337 * @f$\|\boldsymbol{
b}-\boldsymbol{
A}\boldsymbol{c}\|<10^{-8}\|\boldsymbol{
b}\|@f$.
2338 * The maximum number of iteration steps is
set to `system_rhs.size()` as
2339 * the conjugate
gradient algorithm is supposed to find the solution in at
2340 * most @f$m@f$ steps
for an @f$m \times m@f$ system
matrix. This function also
2341 * distributes constraints. The constraints are only used to enforce the
2342 * Dirichlet boundary condition.
2345 *
void BaseSolver::solve()
2347 *
SolverControl control(system_rhs.size(), 1e-8 * system_rhs.l2_norm());
2353 *
preconditioner.
initialize(system_matrix, 1.2);
2355 *
cg.solve(system_matrix, solution, system_rhs, preconditioner);
2357 *
constraints.distribute(solution);
2362 * The following two
functions compute the error norms and
project the exact
2366 *
void BaseSolver::compute_error_norms()
2368 *
if (exact_solution)
2376 *
dof_handler.get_fe(0).degree + 4),
2385 *
void BaseSolver::project_exact_solution_fcn()
2387 *
if (Settings::project_exact_solution && exact_solution)
2391 *
constraints_empty.
clear();
2395 *
constraints_empty.close();
2399 *
constraints_empty,
2401 *
(dof_handler.get_fe().degree + 1) :
2402 *
(dof_handler.get_fe().degree + 2)),
2404 *
projected_exact_solution);
2410 * The following function saves the solution and the @f$L_2@f$ error
norm into a
2411 * `.vtu` file. If `Settings::project_exact_solution =
true`, the projected
2412 * exact solution is saved as well.
2415 *
void BaseSolver::output_results() const
2417 *
const bool fe_is_vector = ((stage == 2) || (stage == 3));
2419 *
const std::string name = (fe_is_vector ?
"VectorField" :
"ScalarField");
2420 *
const unsigned int components = (fe_is_vector ? 2 : 1);
2422 *
component_interpretation =
2427 *
std::vector<std::string> solution_names(components, name);
2428 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
2429 *
data_component_interpretation(components, component_interpretation);
2436 *
data_component_interpretation);
2438 *
if (exact_solution)
2440 *
data_out.add_data_vector(L2_per_cell,
"L2norm");
2442 *
if (Settings::project_exact_solution)
2444 *
std::vector<std::string> solution_names_exact(components,
2446 *
data_out.add_data_vector(dof_handler,
2447 *
projected_exact_solution,
2448 *
solution_names_exact,
2449 *
data_component_interpretation);
2455 *
data_out.set_flags(flags);
2459 *
data_out.build_patches(mapping,
2460 *
dof_handler.get_fe(0).degree + 2,
2463 *
std::ofstream ofs(file_name +
".vtu");
2464 *
data_out.write_vtu(ofs);
2469 * The following function clears the memory
for the next solver.
2472 *
void BaseSolver::clear()
2474 *
system_matrix.clear();
2475 *
system_rhs.reinit(0);
2484 *
void BaseSolver::run()
2490 *
compute_error_norms();
2491 *
project_exact_solution_fcn();
2498 * The following is the straightforward implementation of the get
functions.
2501 *
double BaseSolver::get_L2_norm() const
2506 *
unsigned int BaseSolver::get_n_cells() const
2508 *
return triangulation_rhs.n_active_cells();
2513 *
return dof_handler.n_dofs();
2518 *
return triangulation_rhs;
2523 *
return dof_handler;
2535 * <a name=
"step_98-SolverT"></a>
2536 * <h3>Solver -
T</h3>
2540 * The following
namespace contains the code related to the computation of the
2541 * current vector potential, @f$T@f$.
2546 *
using namespace BaseClasses;
2550 * We derive the solver from the `BaseSolver`
class. What is left to
do is
2551 * to initialize the `BaseSolver`,
override the two
virtual functions
2552 * (`make_mesh` and `system_matrix_local`), and implement the function
2553 * `free_current_density`. The function `free_current_density` implements
2554 * the closed-form analytical expression
for @f$\vec{
J}_f@f$ on the right-hand
2555 * side of the [div-grad equation](@ref Step98_FourStages).
2558 *
class Solver : public BaseSolver
2561 *
Solver() = delete;
2562 *
Solver(
const unsigned int p,
2563 *
const unsigned int r,
2564 *
const unsigned int mapping_degree,
2565 *
const std::string &file_name =
"data",
2569 *
virtual void make_mesh() override final;
2571 *
system_matrix_local(const CellIteratorPair &IP,
2572 *
AssemblyScratchData &scratch_data,
2573 *
AssemblyCopyData ©_data) override final;
2576 *
const
unsigned int refinement_parameter;
2578 *
void free_current_density(const
std::vector<
Point<2>> &p,
2579 *
const
types::material_id material_id,
2580 *
std::vector<
Tensor<1, 2>> &values) const;
2585 * Following is the implementation of the constructor.
2586 * We use the
first constructor of the `BaseSolver` class as
2587 * the solver is used at the
2588 * [
first stage](@ref Step98_FourStages). By looking at the
2589 * [expressions](@ref Step98_Numerical_Recipe_T)
2590 *
for @f$A_{ij}@f$ and @f$b_i@f$ we can conclude that to compute them we need
2592 * weights multiplied by the Jacobian
determinant(`JxW`). The quadrature
2593 * points are needed to sample the closed-form analytical expression
for
2594 * @f$\vec{
J}_f@f$. Accordingly, we use the update flags
2596 *
for the
FE_Q finite elements. This time we
do not use the finite elements
2597 * that model the potential on the right-hand side of the equation, so we
2601 *
Solver::Solver(
const unsigned int p,
2602 *
const unsigned int r,
2603 *
const unsigned int mapping_degree,
2604 *
const std::string &file_name,
2606 *
: BaseSolver(mapping_degree,
2614 *
, refinement_parameter(r)
2619 * The following function loads the mesh, creates manifolds, bounds the
2620 * manifolds to the manifold IDs, refines the mesh, and distributes the dofs.
2621 * Note that we are allowed to create the manifolds locally in the function as
2622 * the triangulation
object keeps copies of the manifolds, see @ref step_65
"step-65". Also
2623 * recall that we have shifted the task of distributing the dofs from the
2624 * `setup` function to the `make_mesh` function to evade the necessity to make
2625 * the `setup` function
virtual. This allows us to keep
one `setup` function
2626 * in the `BaseSolver`
class that serves the needs of all derived classes.
2629 *
void Solver::make_mesh()
2634 *
gridin.read_msh(
"circle.msh");
2642 *
transfinite_manifold.
initialize(triangulation);
2643 *
triangulation.set_manifold(transfinite_interpolation_manifold_id,
2644 *
transfinite_manifold);
2646 *
triangulation.refine_global(refinement_parameter);
2648 *
dof_handler.reinit(triangulation);
2649 *
dof_handler.distribute_dofs(fe);
2654 * The following function assembles a fraction of
2655 * [the system
matrix and the system right-hand side](@ref Step98_Numerical_Recipe_T)
2656 * related to a single cell. These fractions are
2657 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied to
2658 * `system_matrix` and `system_rhs` by
WorkStream.
2661 *
void Solver::system_matrix_local(
const CellIteratorPair &IP,
2662 *
AssemblyScratchData &scratch_data,
2663 *
AssemblyCopyData ©_data)
2667 *
copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
2668 *
scratch_data.dofs_per_cell);
2670 *
copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
2672 *
copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
2676 *
scratch_data.fe_values.reinit(cell);
2678 *
Solver::free_current_density(scratch_data.fe_values.get_quadrature_points(),
2679 *
cell->material_id(),
2680 *
scratch_data.vectors_list_rhs);
2682 *
for (
unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index)
2684 *
for (
unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i)
2686 *
for (
unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j)
2688 *
copy_data.cell_matrix(i, j) +=
2689 *
(scratch_data.fe_values[se].gradient(
2691 *
scratch_data.fe_values[se].gradient(
2694 *
scratch_data.fe_values.JxW(q_index);
2696 *
copy_data.cell_rhs(i) +=
2697 *
(scratch_data.vectors_list_rhs[q_index][0] *
2698 *
scratch_data.fe_values[se].gradient(i, q_index)[1] -
2699 *
scratch_data.vectors_list_rhs[q_index][1] *
2700 *
scratch_data.fe_values[se].gradient(i,
2702 *
scratch_data.fe_values.JxW(
2707 *
cell->get_dof_indices(copy_data.local_dof_indices);
2712 * The following function implements the closed form analytical
2713 * [expression](@ref Step98_Equation_1_Jf)
2714 *
for @f$\vec{
J}_f@f$ on the right-hand side of the
2715 * [div-grad equation](@ref Step98_PDE_T).
2718 *
void Solver::free_current_density(
const std::vector<
Point<2>> &p,
2723 *
ExcDimensionMismatch(p.size(),
values.size()));
2725 *
if ((material_id == Settings::material_id_free_space) ||
2726 *
(material_id == Settings::material_id_core))
2729 *
if (material_id == Settings::material_id_free_current)
2730 *
for (
unsigned int i = 0; i <
values.size(); i++)
2732 *
ExactSolutions::volume_free_current_density(p[i], Settings::K0);
2739 * <a name=
"step_98-ProjectorfromHgradtoHdiv"></a>
2740 * <h3>Projector from H(grad) to H(div)</h3>
2744 * The following
namespace contains all the code related to the computation of
2745 * the free-current density, @f$\vec{
J}_f@f$.
2748 *
namespace ProjectorHgradToHdiv
2750 *
using namespace BaseClasses;
2754 * We derive the solver from the `BaseSolver`
class. What is left to
do is
2755 * to initialize the `BaseSolver` and
override the two
virtual functions
2756 * (`make_mesh` and `system_matrix_local`).
2759 *
class Solver : public BaseSolver
2762 *
Solver() = delete;
2763 *
Solver(
const unsigned int p,
2765 *
const unsigned int mapping_degree,
2769 *
const std::string &file_name =
"data",
2773 *
virtual void make_mesh() override final;
2775 *
system_matrix_local(const CellIteratorPair &IP,
2776 *
AssemblyScratchData &scratch_data,
2777 *
AssemblyCopyData ©_data) override final;
2784 * Following is the implementation of the constructor.
2785 * We use the
second constructor of the `BaseSolver` class as
2786 * the solver is used at the
2787 * [
second stage](@ref Step98_FourStages). By looking at the
2788 * [expressions](@ref Step98_Numerical_Recipe_Jf)
2789 *
for @f$A_{ij}@f$ and @f$b_i@f$ we can conclude that to compute them we need
2790 *
values of the shape
functions and the quadrature weights multiplied by the
2794 * there is a numerically computed potential, @f$T@f$, on the right-hand side of
2795 * the [equation](@ref Step98_PDE_Jf). It is modeled by the
FE_Q finite elements.
2796 * To compute the right-hand side, we need gradients of the shape
functions.
2801 *
Solver::Solver(
const unsigned int p,
2802 *
const unsigned int mapping_degree,
2806 *
const std::string &file_name,
2808 *
: BaseSolver(triangulation_rhs,
2821 * At the
second stage we
do not
load the mesh. We reuse the mesh loaded
2822 * at the
first stage. Consequently, we just need to distribute the dofs.
2825 *
void Solver::make_mesh()
2827 *
dof_handler.reinit(triangulation_rhs);
2828 *
dof_handler.distribute_dofs(fe);
2833 * The following function assembles a fraction of
2834 * [the system
matrix and the system right-hand side](@ref Step98_Numerical_Recipe_Jf)
2835 * related to a single cell. These fractions are
2836 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied to
2837 * `system_matrix` and `system_rhs` by
WorkStream.
2840 *
void Solver::system_matrix_local(
const CellIteratorPair &IP,
2841 *
AssemblyScratchData &scratch_data,
2842 *
AssemblyCopyData ©_data)
2846 *
copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
2847 *
scratch_data.dofs_per_cell);
2849 *
copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
2851 *
copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
2857 *
scratch_data.fe_values.reinit(cell);
2858 *
scratch_data.fe_values_rhs.reinit(cell_rhs);
2860 *
scratch_data.fe_values_rhs.get_function_gradients(
2861 *
scratch_data.dofs_rhs, scratch_data.vectors_list_rhs);
2863 *
for (
unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index)
2865 *
for (
unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i)
2867 *
for (
unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j)
2869 *
copy_data.cell_matrix(i, j) +=
2870 *
scratch_data.fe_values[ve].value(i, q_index) *
2871 *
scratch_data.fe_values[ve].value(j,
2873 *
scratch_data.fe_values.JxW(q_index);
2876 *
copy_data.cell_rhs(i) +=
2877 *
(scratch_data.vectors_list_rhs[q_index][1] *
2878 *
scratch_data.fe_values[ve].value(i, q_index)[0] -
2879 *
scratch_data.vectors_list_rhs[q_index][0] *
2880 *
scratch_data.fe_values[ve].value(i, q_index)[1]) *
2881 *
scratch_data.fe_values.JxW(
2886 *
cell->get_dof_indices(copy_data.local_dof_indices);
2895 * <a name=
"step_98-SolverA"></a>
2896 * <h3>Solver -
A</h3>
2900 * The following
namespace contains all the code related to the computation of
2901 * the magnetic vector potential, @f$\vec{
A}@f$.
2906 *
using namespace BaseClasses;
2910 * We derive the solver from the `BaseSolver`
class. What is left to
do is
2911 * to initialize the `BaseSolver`,
override the two
virtual functions
2912 * (`make_mesh` and `system_matrix_local`), and implement the function
2916 *
class Solver : public BaseSolver
2919 *
Solver() = delete;
2920 *
Solver(
const unsigned int p,
2921 *
const unsigned int mapping_degree,
2925 *
const std::string &file_name =
"data",
2929 *
virtual void make_mesh() override final;
2931 *
system_matrix_local(const CellIteratorPair &IP,
2932 *
AssemblyScratchData &scratch_data,
2933 *
AssemblyCopyData ©_data) override final;
2937 *
void permeability(const
types::material_id material_id,
2938 *
std::vector<
double> &values) const;
2943 * Following is the implementation of the constructor.
2944 * We use the
second constructor of the `BaseSolver` class as
2945 * the solver is used at the
2946 * [third stage](@ref Step98_FourStages). By looking at the
2947 * [expressions](@ref Step98_Numerical_Recipe_A)
2948 *
for @f$A_{ij}@f$ and @f$b_i@f$ we can conclude that to compute them we need
2951 * finite elements. Accordingly, we use the update flags `
update_values`,
2953 * elements. This time there is a numerically computed potential, @f$T@f$, on
2954 * the right-hand side of the
2955 * [equation](@ref Step98_PDE_A). It is modeled by the
FE_Q finite elements.
2956 * To compute it, we need values of the shape
functions. Accordingly, we use
2960 *
Solver::Solver(
const unsigned int p,
2961 *
const unsigned int mapping_degree,
2965 *
const std::string &file_name,
2967 *
: BaseSolver(triangulation_rhs,
2981 * At the third stage we
do not
load the mesh. We reuse the mesh loaded at
2982 * the
first stage. Consequently, we just need to distribute the dofs.
2985 *
void Solver::make_mesh()
2987 *
dof_handler.reinit(triangulation_rhs);
2988 *
dof_handler.distribute_dofs(fe);
2993 * The following function assembles a fraction of
2994 * [the system
matrix and the system right-hand side](@ref Step98_Numerical_Recipe_A)
2995 * related to a single cell. These fractions are
2996 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied to
2997 * `system_matrix` and `system_rhs` by
WorkStream.
3000 *
void Solver::system_matrix_local(
const CellIteratorPair &IP,
3001 *
AssemblyScratchData &scratch_data,
3002 *
AssemblyCopyData ©_data)
3006 *
copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
3007 *
scratch_data.dofs_per_cell);
3009 *
copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
3011 *
copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
3017 *
scratch_data.fe_values.reinit(cell);
3018 *
scratch_data.fe_values_rhs.reinit(cell_rhs);
3020 *
Solver::permeability(cell->material_id(), scratch_data.permeability_list);
3022 *
scratch_data.fe_values_rhs.get_function_values(
3023 *
scratch_data.dofs_rhs, scratch_data.values_list_rhs);
3025 *
for (
unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index)
3027 *
for (
unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i)
3029 *
for (
unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j)
3031 *
copy_data.cell_matrix(i, j) +=
3032 *
(1.0 / scratch_data.permeability_list[q_index]) *
3033 *
(scratch_data.fe_values[ve].curl(
3035 *
scratch_data.fe_values[ve].curl(
3038 *
Settings::eta_squared *
3039 *
scratch_data.fe_values[ve].value(i,
3041 *
scratch_data.fe_values[ve].value(j, q_index)
3043 *
scratch_data.fe_values.JxW(q_index);
3045 *
copy_data.cell_rhs(i) +=
3046 *
(scratch_data.values_list_rhs[q_index] *
3047 *
scratch_data.fe_values[ve].curl(i, q_index)) *
3048 *
scratch_data.fe_values.JxW(q_index);
3052 *
cell->get_dof_indices(copy_data.local_dof_indices);
3057 * The following function implements the
3058 * [equation](@ref Step98_Equation_MU)
3063 *
std::vector<double> &values)
const
3065 *
if (material_id == Settings::material_id_core)
3066 *
std::fill(
values.begin(),
values.end(), Settings::mu_1);
3068 *
std::fill(
values.begin(),
values.end(), Settings::mu_0);
3075 * <a name=
"step_98-ProjectorfromHcurltoL2"></a>
3076 * <h3>Projector from H(curl) to
L2</h3>
3080 * The following
namespace contains all the code related to the computation of
3081 * the magnetic field, @f$B@f$.
3084 *
namespace ProjectorHcurlToL2
3086 *
using namespace BaseClasses;
3090 * We derive the solver from the `BaseSolver`
class. What is left to
do is
3091 * to initialize the `BaseSolver` and
override the two
virtual functions
3092 * (`make_mesh` and `system_matrix_local`).
3095 *
class Solver : public BaseSolver
3098 *
Solver() = delete;
3099 *
Solver(
const unsigned int p,
3100 *
const unsigned int mapping_degree,
3104 *
const std::string &file_name =
"data",
3108 *
virtual void make_mesh() override final;
3110 *
system_matrix_local(const CellIteratorPair &IP,
3111 *
AssemblyScratchData &scratch_data,
3112 *
AssemblyCopyData ©_data) override final;
3119 * Following is the implementation of the constructor.
3120 * We use the
second constructor of the `BaseSolver` class as
3121 * the solver is used at the
3122 * [fourth stage](@ref Step98_FourStages). By looking at the
3123 * [expressions](@ref Step98_Numerical_Recipe_B)
3124 *
for @f$A_{ij}@f$ and @f$b_i@f$ we can conclude that to compute them we need
3125 *
values of the shape
functions and the quadrature weights multiplied by the
3129 * numerically computed potential, @f$\vec{
A}@f$, on the right-hand side of the
3130 * [equation](@ref Step98_PDE_B). It is modeled by the
FE_Nedelec finite
3131 * elements. To compute the right-hand side, we need gradients of the shape
3136 *
Solver::Solver(
const unsigned int p,
3137 *
const unsigned int mapping_degree,
3141 *
const std::string &file_name,
3143 *
: BaseSolver(triangulation_ext,
3156 * At the fourth stage we
do not
load the mesh. We reuse the mesh loaded
3157 * at the
first stage. Consequently, we just need to distribute the dofs.
3160 *
void Solver::make_mesh()
3162 *
dof_handler.reinit(triangulation_rhs);
3163 *
dof_handler.distribute_dofs(fe);
3168 * The following function assembles a fraction of
3169 * [the system
matrix and the system right-hand side](@ref Step98_Numerical_Recipe_B)
3170 * related to a single cell. These fractions are
3171 * `copy_data.cell_matrix` and `copy_data.cell_rhs`. They are copied to
3172 * `system_matrix` and `system_rhs` by
WorkStream.
3175 *
void Solver::system_matrix_local(
const CellIteratorPair &IP,
3176 *
AssemblyScratchData &scratch_data,
3177 *
AssemblyCopyData ©_data)
3181 *
copy_data.cell_matrix.reinit(scratch_data.dofs_per_cell,
3183 *
scratch_data.dofs_per_cell);
3185 *
copy_data.cell_rhs.reinit(scratch_data.dofs_per_cell);
3187 *
copy_data.local_dof_indices.resize(scratch_data.dofs_per_cell);
3193 *
scratch_data.fe_values.reinit(cell);
3194 *
scratch_data.fe_values_rhs.reinit(cell_rhs);
3196 *
scratch_data.fe_values_rhs.get_function_gradients(
3197 *
scratch_data.dofs_rhs, scratch_data.vectors_vectors_list_rhs);
3199 *
for (
unsigned int q_index = 0; q_index < scratch_data.n_q_points; ++q_index)
3201 *
for (
unsigned int i = 0; i < scratch_data.dofs_per_cell; ++i)
3203 *
for (
unsigned int j = 0; j < scratch_data.dofs_per_cell; ++j)
3205 *
copy_data.cell_matrix(i, j) +=
3206 *
scratch_data.fe_values[se].value(i, q_index) *
3207 *
scratch_data.fe_values[se].value(j, q_index) *
3208 *
scratch_data.fe_values.JxW(q_index);
3211 *
copy_data.cell_rhs(i) +=
3212 *
(scratch_data.vectors_vectors_list_rhs[q_index][1][0] -
3213 *
scratch_data.vectors_vectors_list_rhs[q_index][0][1]) *
3214 *
scratch_data.fe_values[se].value(i, q_index) *
3215 *
scratch_data.fe_values.JxW(
3220 *
cell->get_dof_indices(copy_data.local_dof_indices);
3227 * <a name=
"step_98-Themainloop"></a>
3233 * function. The implementation of the
loop is straightforward -
3234 * we create and
run the four solvers
one-by-
one and
copy the relevant
3235 *
data into convergence tables.
3238 *
class MagneticProblem
3243 *
if (Settings::n_threads_max != 0)
3246 *
MainOutputTable table_T(2);
3247 *
MainOutputTable table_Jf(2);
3248 *
MainOutputTable table_A(2);
3249 *
MainOutputTable table_B(2);
3251 *
std::cout <<
"Solving for (p' = " << Settings::fe_degree + 1
3252 *
<<
"; p = " << Settings::fe_degree <<
"): " << std::flush;
3254 *
for (
unsigned int r = 1; r < 5; r++)
3256 *
table_T.add_value(
"r", r);
3257 *
table_T.add_value(
"p", Settings::fe_degree + 1);
3259 *
table_Jf.add_value(
"r", r);
3260 *
table_Jf.add_value(
"p", Settings::fe_degree);
3262 *
table_A.add_value(
"r", r);
3263 *
table_A.add_value(
"p", Settings::fe_degree);
3265 *
table_B.add_value(
"r", r);
3266 *
table_B.add_value(
"p", Settings::fe_degree);
3270 * Stage 1. Computing @f$T@f$.
3276 *
std::cout <<
"T " << std::flush;
3278 *
ExactSolutions::CurrentVectorPotential T_exact;
3280 *
SolverT::Solver
T(Settings::fe_degree + 1,
3282 *
Settings::mapping_degree,
3283 *
"T_p" + std::to_string(Settings::fe_degree + 1) +
3284 *
"_r" + std::to_string(r),
3289 *
table_T.add_value(
"ndofs",
T.get_n_dofs());
3290 *
table_T.add_value(
"ncells",
T.get_n_cells());
3291 *
table_T.add_value(
"L2",
T.get_L2_norm());
3295 * Stage 2. Computing @f$\vec{
J}_f@f$.
3301 *
std::cout <<
"Jf " << std::flush;
3303 *
ExactSolutions::FreeCurrentDensity Jf_exact;
3305 *
ProjectorHgradToHdiv::Solver Jf(Settings::fe_degree,
3306 *
Settings::mapping_degree,
3308 *
T.get_dof_handler(),
3311 *
std::to_string(Settings::fe_degree) +
3312 *
"_r" + std::to_string(r),
3317 *
table_Jf.add_value(
"ndofs", Jf.get_n_dofs());
3318 *
table_Jf.add_value(
"ncells", Jf.get_n_cells());
3319 *
table_Jf.add_value(
"L2", Jf.get_L2_norm());
3323 * Stage 3. Computing @f$\vec{
A}@f$.
3329 *
std::cout <<
"A " << std::flush;
3331 *
ExactSolutions::MagneticVectorPotential A_exact;
3333 *
SolverA::Solver
A(Settings::fe_degree,
3334 *
Settings::mapping_degree,
3336 *
T.get_dof_handler(),
3338 *
"A_p" + std::to_string(Settings::fe_degree) +
"_r" +
3339 *
std::to_string(r),
3344 *
table_A.add_value(
"ndofs",
A.get_n_dofs());
3345 *
table_A.add_value(
"ncells",
A.get_n_cells());
3346 *
table_A.add_value(
"L2",
A.get_L2_norm());
3350 * Stage 4. Computing @f$B@f$.
3356 *
std::cout <<
"B " << std::flush;
3358 *
ExactSolutions::MagneticField B_exact;
3360 *
ProjectorHcurlToL2::Solver B(Settings::fe_degree,
3361 *
Settings::mapping_degree,
3363 *
A.get_dof_handler(),
3366 *
std::to_string(Settings::fe_degree) +
3367 *
"_r" + std::to_string(r),
3371 *
table_B.add_value(
"ndofs", B.get_n_dofs());
3372 *
table_B.add_value(
"ncells", B.get_n_cells());
3373 *
table_B.add_value(
"L2", B.get_L2_norm());
3381 *
table_T.save(
"table_T_p" + std::to_string(Settings::fe_degree + 1));
3382 *
table_Jf.save(
"table_Jf_p" + std::to_string(Settings::fe_degree));
3383 *
table_A.save(
"table_A_p" + std::to_string(Settings::fe_degree));
3384 *
table_B.save(
"table_B_p" + std::to_string(Settings::fe_degree));
3385 *
std::cout << std::endl;
3389 *
int main(
int argc,
char **argv)
3395 * gmsh has global state, so
set that up in the normal way:
3401 *
MagneticProblem problem;
3404 *
catch (std::exception &exc)
3406 *
std::cerr << std::endl
3408 *
<<
"----------------------------------------------------"
3410 *
std::cerr <<
"Exception on processing: " << std::endl
3411 *
<< exc.what() << std::endl
3412 *
<<
"Aborting!" << std::endl
3413 *
<<
"----------------------------------------------------"
3420 *
std::cerr << std::endl
3422 *
<<
"----------------------------------------------------"
3424 *
std::cerr <<
"Unknown exception!" << std::endl
3425 *
<<
"Aborting!" << std::endl
3426 *
<<
"----------------------------------------------------"
3434<a name=
"step_98-Results"></a><h1>Results</h1>
3437The program generates the following output in the command line
interface by
3441Solving
for (p' = 1; p = 0):
T Jf
A B
T Jf
A B
T Jf
A B
T Jf
A B
3444The program assumes the finite elements of the lowermost degree, i.
e., @f$p
' = 1@f$
3445for the FE_Q finite elements and @f$p=0@f$ for other finite elements. To change the
3446degree of the finite elements, say @f$p' = 3@f$ and @f$p = 2@f$,
one needs to change
3447the setting `
Settings::fe_degree = 2` and rebuild the program. The degree of the
3448FE_Q finite elements will be computed automatically as @f$p
'= p + 1@f$.
3450The program dumps a number of files in the current directory. In the default
3451configuration these files are:
3452- `.vtu` files. They contain the computed vector fields. Recall that the spherical
3453 manifold and transfinite interpolation manifold are attached to many cell faces.
3454 Consequently, these cell faces are curved. Furthermore, the shape functions are
3455 mapped from the reference cell to the real mesh cells by the second-order mapping
3456 to accommodate the cells with curved faces. For these reasons, one needs to use a
3457 visualization software that can deal with curved faces and the higher-order mapping.
3458 A fresh version of ParaView is recommended. Visit did not have this feature
3459 at the time this tutorial was written (in early 2026). The
3460 <a href="https://github.com/dealii/dealii/wiki/Notes-on-visualizing-high-order-output"> Notes on visualizing high order output</a>
3461 provide more information on this topic.
3462- `.tex` files. These files contain the convergence tables.
3464The following provides examples of the convergence tables simulated with the
3465default settings for three different degrees of the finite elements,
3466@f$p = 0, 1, 2@f$ (recall that @f$p' = p + 1@f$).
3469<caption>Convergence table @f$T@f$.</caption>
3475 <th>@f$\|e\|_{L^2}@f$</th>
3476 <th>@f$\alpha_{L^2}@f$</th>
3577<caption>Convergence table @f$\vec{J}_f@f$.</caption>
3583 <th>@f$\|e\|_{L^2}@f$</th>
3584 <th>@f$\alpha_{L^2}@f$</th>
3685<caption> Convergence table @f$B@f$.</caption>
3691 <th>@f$\|e\|_{L^2}@f$</th>
3692 <th>@f$\alpha_{L^2}@f$</th>
3792The following notations were used in the headers of the tables:
3794- p', p - the degree of the finite elements.
3796- r - the mesh refinement parameter, i.
e., the number of global mesh
3799- cells - the total number of active cells.
3801- dofs - the number of degrees of freedom.
3803-@f$\|
e\|_{
L^2}@f$ - the @f$L^2@f$ error
norm.
3805-@f$\alpha_{
L^2}@f$ - the order of convergence of the @f$L^2@f$ error
norm.
3807Let us contemplate these convergence tables
for a brief moment. The
first table
3808illustrates convergence of the numerically computed current vector potential,
3809@f$T@f$. In this particular case we can expect the order of the convergence rate of
3810@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
3811elements of the lowermost degree are used, @f$p'=1@f$, the order of
3812convergence rate is at the upper boundary of the expected
values,
3813@f$\alpha_{
L^2} \approx 2.0@f$. If the finite elements of the
second degree are
3814used, @f$p
'=2@f$, the order of the convergence rate is higher than expected,
3815@f$\alpha_{L^2} \approx 4.0@f$. This means that the @f$L^2@f$ error norm converges to
3816zero at a rate higher than theoretically possible. Most likely this is due to
3817the fact that the current vector potential has a particularly simple form,
3818[Equation 3](@ref Step98_Equation_3_T). It is either constant or changes as
3819the second-order monomial,
3823The second order polynomial approximates this behavior exactly. That is to
3824say, in this particular case we have a lucky situation in which the shape
3825functions can approximate the numerical solution exactly within each mesh
3826cell. This, most likely, explains the extra rapid convergence of the @f$L^2@f$
3827norm. Note also that the error norms, @f$\|e\|_{L^2}@f$, at @f$p'=2@f$ and @f$p
'=3@f$
3828are the same for the same values of the mesh refinement parameter, @f$r@f$. This
3829is, most likely, due to the fact that the second-order shape functions,
3830@f$p'=2@f$,
approximate the solution exactly within each mesh cell and the
3831third-order monomials of the shape function at @f$p
'=3@f$ have absolutely
3832nothing to contribute to the quality of approximation.
3834The second table illustrates convergence of the numerically computed
3835free-current density, @f$\vec{J}_f@f$. The free-current density has been
3836computed as a derivative of current vector potential,
3838\vec{J}_f = \vec{\nabla}\overset{V} {\times} T.
3840The derivative reduces the order of the convergence rate by one. Most
3841likely in this particular case the error made in computing @f$\vec{J}_f@f$
3842is defined by the error made in computing @f$T@f$. For this reason, the order
3843of convergence in the second table equals the order of convergence in
3844the first table minus one.
3846Due to the implicit gauge we cannot observe the convergence of the magnetic
3847vector potential, @f$\vec{A}@f$. Instead, we can observe the convergence on the
3848magnetic field, @f$B@f$, given in the third table. It follows the same pattern:
3849The rate of convergence at the lower degrees of the finite elements is at
3850the best expected value, @f$\alpha_{L^2} = p + 1@f$; the rate of convergence at
3851the higher degrees of the finite elements is better than expected. The extra
3852rapid convergence at the higher degrees of the finite elements is explained
3853by the relatively simple form of the field being approximated.
3855The figures below illustrate the current vector potential, @f$T@f$, the free-current
3856density, @f$\vec{J}_f@f$, and magnetic field, @f$B@f$, computed with the following
3857settings: @f$p = 2@f$ and @f$r = 4@f$. Visual inspection of the magnetic vector
3858potential, @f$\vec{A}@f$, is not very informative as its conservative portion is
3863 <img src="https://dealii.org/images/steps/developer/step-98-result-T.svg"
3864 alt="The result - current vector potential" height="531">
3870 <img src="https://dealii.org/images/steps/developer/step-98-result-Jf.svg"
3871 alt="The result - free current density" height="531">
3877 <img src="https://dealii.org/images/steps/developer/step-98-result-B.svg"
3878 alt="The result - magnetic field" height="531">
3882The images above suggest that the computed fields do not exhibit any irregular
3883behavior (the first image in the table below illustrates how irregular
3884behavior can look like). The fields on these images closely resemble the
3885corresponding closed-form analytical expressions given in the introduction.
3886From the first glans the convergence tables above may appear somewhat strange
3887due to extra rapid convergence at the higher degrees of the finite elements.
3888One, however, can argue that the extra rapid convergence can be explained by
3889the simple polynomial form of the fields being approximated. One thing is
3890certain - the convergence rates presented in these tables are at the best
3891theoretically expected values or better.
3893@anchor Step98_PossibilitiesForExtensions
3894<a name="step_98-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
3896Let us consider the two-dimensional curl-curl partial differential equation
3899\vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu} \vec{\nabla}\overset{S}{\times}\vec{A}\bigg)
3900+ \eta^2 \vec{A} = \vec{\nabla}\overset{V}{\times} T.
3902We can compute the divergence of this expression and rearrange the terms as the
3905\vec{\nabla}\cdot \vec{A} = \frac{1}{\eta^2} \vec{\nabla}\cdot
3906\bigg[\vec{\nabla}\overset{V}{\times} T -
3907\vec{\nabla}\overset{V}{\times}\bigg(\dfrac{1}{\mu}
3908\vec{\nabla}\overset{S}{\times}\vec{A}\bigg)\bigg].
3910The divergence of the vector curl equals zero, see the introduction. Therefore,
3911the right-hand side of the last equation equals zero,
3913\vec{\nabla}\cdot \vec{A} = 0.
3915That is to say, the @f$\eta^2@f$ gauging term can be considered to be the Coulomb
3916gauge at least in theory. In practice the situation is a bit more complicated.
3917What the @f$\eta^2@f$ gauging term does depends on the value of @f$\eta^2@f$. The table
3918below attempts to express this very point.
3921 <img src="https://dealii.org/images/steps/developer/step-98-eta.svg"
3922 alt="The result - eta-squared table" height="1062">
3926This table presents four simulations for four different values of @f$\eta^2@f$. In
3927all four simulations the degree of the finite elements and the mesh refinement
3928parameter were @f$p=2@f$ and @f$r=1@f$, respectively. The setting
3929`Settings::project_exact_solution` was set to `true`.
3931At the setting @f$\eta^2=\dfrac{10^{-10}}{\mu_0}@f$ the computed magnetic vector
3932potential (the orange curve on the @f$|\vec{A}(x,0)|@f$ plot) looks exactly the same
3933as the projected exact solution gauged by the Coulomb gauge (the blue curve). At
3934this value of @f$\eta^2@f$ the term @f$\eta^2\vec{A}@f$ acts like the Coulomb gauge.
3935At the default setting, @f$\eta^2=0@f$, the solution, @f$\vec{A}@f$, is contaminated by
3936an unknown conservative vector field. The computed magnetic field, @f$B@f$, at this
3937setting looks exactly as the corresponding exact solution because the
3938conservative portion of @f$\vec{A}@f$ is filtered out by the process of computing
3939the magnetic field, @f$B = \vec{\nabla}\overset{S}{\times}\vec{A}@f$. At the setting
3940@f$\eta^2=\dfrac{10^{-6}}{\mu_0}@f$ the solenoidal part of the computed magnetic
3941vector potential deviates from the exact solution. We can deduce this from the
3942fact that the computed magnetic
3943field deviates from the exact expression for the magnetic field on the @f$B(x,0)@f$
3944plot. The error in the solenoidal part of the computed @f$\vec{A}@f$ is due to
3945the fact that introduction of the gauging term, @f$\eta^2\vec{A}@f$, modifies the
3947curl-curl equation. So, strictly speaking, we are solving a different partial
3948differential equation. If @f$\eta^2@f$ is small, the difference between the initial
3949curl-curl equation and the curl-curl equation modified by adding the gauging
3950term is negligible. Consequently, the error in the solenoidal part of the
3951computed @f$\vec{A}@f$ is negligible as well. Evidently, in this particular case
3952"small" means @f$\eta^2 \ll \dfrac{10^{-6}}{\mu_0}@f$.
3954Normally, we are interested in measurable fields such as magnetic field and
3955consider the magnetic vector potential as a useful tool for computing measurable
3956fields. We can tolerate the presence of an unknown conservative component in the
3957magnetic vector potential, i.e., the situation illustrated by the first two rows
3958in the table above. In such disposition we need to keep @f$\eta^2@f$ as small as
3959possible, i.e., as far away as possible from the situation shown in last row of
3960the table. In this tutorial program we set @f$\eta^2@f$ to zero and increase it just
3961a bit in if the conjugate gradient algorithm cannot converge.
3963Suppose for a moment that we would like to have the solution to the curl-curl
3964equation in terms of a purely solenoidal magnetic vector potential, @f$\vec{A}@f$,
3965that is, a solution gauged by the Coulomb gauge. To get such a solution we need
3966to tweak the @f$\eta^2@f$ parameter. By contemplating
3967the table above one can hypothesize that there exists an optimal value of the
3968gauging parameter, @f$\eta^2_\text{opt}@f$, at which the @f$L^2@f$ error norm computed for
3969@f$\vec{A}@f$ is minimal. The optimal value should be somewhere in between
3970@f$\eta^2=\dfrac{10^{-12}}{\mu_0}@f$ and @f$\eta^2=\dfrac{10^{-6}}{\mu_0}@f$, based on the
3971experiments above. Try to
3972verify this hypothesis by finding the exact value of @f$\eta^2_\text{opt}@f$.
3974In this instance we have a close-form analytical expression of the exact
3975solution, i.e., the expression of @f$\vec{A}@f$ given in the introduction. In a
3976real-life simulation there is no expression of the exact solution. Try to
3977think of a method of blind (meaning without the exact solution) estimation of
3978@f$\eta^2_\text{opt}@f$. Try to implement and test your ideas.
3980Adding the gauging term, @f$\eta^2\vec{A}@f$, to the curl-curl equation converts
3981a positive semidefinite system matrix into a positive definite matrix. Strictly
3982speaking, adding the gauging term modifies the initial curl-curl equation. For
3983this reason, it is important to keep @f$\eta^2@f$ small so the solution is not
3984afflicted by the error induced by adding the gauging term. The definition of
3985"small" here is a bit fuzzy. In absence of the exact solution setting @f$\eta^2@f$
3986at the acceptable level of the error in the solenoidal component of @f$\vec{A}@f$
3987is difficult. It is better to discard the gauging term and use a more
3988sophisticated linear solver. The hypre AMS @cite hypre1998b can solve the
3989systems of linear equations yielded by the curl-curl equation without the
3990@f$\eta^2@f$ gauging term, (@f$\beta=0@f$ in @cite hypre1998b). Try to implement the
3994<a name="step_98-PlainProg"></a>
3995<h1> The plain program</h1>
3996@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