I ended up getting it working, so in case anyone else is interested in
doing something similar, I have attached updated versions of the previous
scripts. I further modified set_periodicity_constraints so that it can take
a function to set the inhomogeneity based on the values of the two support
points associated with the two degrees of freedom that are being tied
together, although it would need to be modified for more general use (the
functions works on the difference of the two points, and further it uses
only the absolute value of the difference in the x-component). I have also
not made changes to allow for adaptive mesh refinement. I think this is a
pretty ugly solution to the general problem of inhomogeneous periodic
boundary constraints, but I couldn't see how to do it without modifying
these functions without largely replicating them.
The issues mentioned in the previous post regarding
distribute_local_to_global() and AffineConstraints::distribute() stemmed
from a misunderstanding of exactly what these functions were doing, and
have been fixed. It does seem that with these constraints, SolverCG does
not work as a solver (I suppose the matrix is not symmetric anymore?), but
either SolverGMRES or the direct solver SparseDirectUMFPACK work.
The script has also been improved based on other tutorials, and now
includes parsing for a parameter file. I tried to include an example
parameter file but I guess the mailing list has a restriction on allowable
file types.
Best,
Alex
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/0c692129-a363-4d67-98a7-8c7e04ba3892n%40googlegroups.com.
#ifndef modded_periodic_functions_h
#define modded_periodic_functions_h
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/utilities.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/fe/mapping_q1.h>
using namespace dealii;
template <typename FaceIterator, typename number, int dim>
void set_periodicity_constraints(
const FaceIterator& face_1,
const typename identity<FaceIterator>::type& face_2,
const FullMatrix<double>& transformation,
AffineConstraints<number>& affine_constraints,
const std::vector<Point<dim>> dofs_to_supports,
Function<dim>& boundary_function,
const ComponentMask& component_mask,
const bool face_orientation,
const bool face_flip,
const bool face_rotation,
const number periodicity_factor) {
static const int spacedim = FaceIterator::AccessorType::space_dimension;
// we should be in the case where face_1 is active, i.e. has no children:
Assert(!face_1->has_children(), ExcInternalError());
Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
// TODO: the implementation makes the assumption that all faces have the
// same number of dofs
// AssertDimension(
// face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(),
// 1);
// AssertDimension(
// face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(),
// 1);
// const unsigned int face_no = 0;
// If face_2 does have children, then we need to iterate over these
// children and set periodic constraints in the inverse direction:
if (face_2->has_children()) {
Assert(face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
ExcNotImplemented());
const unsigned int dofs_per_face =
face_1->get_fe(face_1->nth_active_fe_index(0))
.n_dofs_per_face();
FullMatrix<double> child_transformation(dofs_per_face, dofs_per_face);
FullMatrix<double> subface_interpolation(dofs_per_face, dofs_per_face);
for (unsigned int c = 0; c < face_2->n_children(); ++c) {
// get the interpolation matrix recursively from the one that
// interpolated from face_1 to face_2 by multiplying from the left
// with the one that interpolates from face_2 to its child
const auto& fe = face_1->get_fe(face_1->nth_active_fe_index(0));
fe.get_subface_interpolation_matrix(fe, c, subface_interpolation);
subface_interpolation.mmult(child_transformation, transformation);
set_periodicity_constraints(
face_1,
face_2->child(c),
child_transformation,
affine_constraints,
dofs_to_supports,
boundary_function,
component_mask,
face_orientation,
face_flip,
face_rotation,
periodicity_factor);
}
return;
}
// If we reached this point then both faces are active. Now all
// that is left is to match the corresponding DoFs of both faces.
const unsigned int face_1_index = face_1->nth_active_fe_index(0);
const unsigned int face_2_index = face_2->nth_active_fe_index(0);
Assert(face_1->get_fe(face_1_index) == face_2->get_fe(face_2_index),
ExcMessage("Matching periodic cells need to use the same finite "
"element"));
const FiniteElement<dim, spacedim>& fe = face_1->get_fe(face_1_index);
Assert(component_mask.represents_n_components(fe.n_components()),
ExcMessage("The number of components in the mask has to be either "
"zero or equal to the number of components in the finite "
"element."));
const unsigned int dofs_per_face = fe.n_dofs_per_face();
std::vector<types::global_dof_index> dofs_1(dofs_per_face);
std::vector<types::global_dof_index> dofs_2(dofs_per_face);
face_1->get_dof_indices(dofs_1, face_1_index);
face_2->get_dof_indices(dofs_2, face_2_index);
// If either of the two faces has an invalid dof index, stop. This is
// so that there is no attempt to match artificial cells of parallel
// distributed triangulations.
//
// While it seems like we ought to be able to avoid even calling
// set_periodicity_constraints for artificial faces, this situation
// can arise when a face that is being made periodic is only
// partially touched by the local subdomain.
// make_periodicity_constraints will be called recursively even for
// the section of the face that is not touched by the local
// subdomain.
//
// Until there is a better way to determine if the cells that
// neighbor a face are artificial, we simply test to see if the face
// does not have a valid dof initialization.
for (unsigned int i = 0; i < dofs_per_face; i++)
if (dofs_1[i] == numbers::invalid_dof_index ||
dofs_2[i] == numbers::invalid_dof_index) {
return;
}
// Well, this is a hack:
//
// There is no
// face_to_face_index(face_index,
// face_orientation,
// face_flip,
// face_rotation)
// function in FiniteElementData, so we have to use
// face_to_cell_index(face_index, face
// face_orientation,
// face_flip,
// face_rotation)
// But this will give us an index on a cell - something we cannot work
// with directly. But luckily we can match them back :-]
std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
// Build up a cell to face index for face_2:
for (unsigned int i = 0; i < dofs_per_face; ++i) {
const unsigned int cell_index = fe.face_to_cell_index(
i,
0, /* It doesn't really matter, just
* assume we're on the first face...
*/
true,
false,
false // default orientation
);
cell_to_rotated_face_index[cell_index] = i;
}
// Loop over all dofs on face 2 and constrain them against all
// matching dofs on face 1:
for (unsigned int i = 0; i < dofs_per_face; ++i) {
auto component_index {fe.face_system_to_component_index(i).first};
// Obey the component mask
if ((component_mask.n_selected_components(fe.n_components()) !=
fe.n_components()) &&
!component_mask[component_index])
continue;
// We have to be careful to treat so called "identity
// constraints" special. These are constraints of the form
// x1 == constraint_factor * x_2. In this case, if the constraint
// x2 == 1./constraint_factor * x1 already exists we are in trouble.
//
// Consequently, we have to check that we have indeed such an
// "identity constraint". We do this by looping over all entries
// of the row of the transformation matrix and check whether we
// find exactly one nonzero entry. If this is the case, set
// "is_identity_constrained" to true and record the corresponding
// index and constraint_factor.
bool is_identity_constrained = false;
unsigned int target = numbers::invalid_unsigned_int;
number constraint_factor = periodicity_factor;
constexpr double eps = 1.e-13;
for (unsigned int jj = 0; jj < dofs_per_face; ++jj) {
const auto entry = transformation(i, jj);
if (std::abs(entry) > eps) {
if (is_identity_constrained) {
// We did encounter more than one nonzero entry, so
// the dof is not identity constrained. Set the
// boolean to false and break out of the for loop.
is_identity_constrained = false;
break;
}
is_identity_constrained = true;
target = jj;
constraint_factor = entry * periodicity_factor;
}
}
// Next, we work on all constraints that are not identity
// constraints, i.e., constraints that involve an interpolation
// step that constrains the current dof (on face 2) to more than
// one dof on face 1.
if (!is_identity_constrained) {
// The current dof is already constrained. There is nothing
// left to do.
if (affine_constraints.is_constrained(dofs_2[i]))
continue;
// Enter the constraint piece by piece:
affine_constraints.add_line(dofs_2[i]);
for (unsigned int jj = 0; jj < dofs_per_face; ++jj) {
// Get the correct dof index on face_1 respecting the
// given orientation:
const unsigned int j =
cell_to_rotated_face_index[fe.face_to_cell_index(
jj,
0,
face_orientation,
face_flip,
face_rotation)];
if (std::abs(transformation(i, jj)) > eps)
affine_constraints.add_entry(
dofs_2[i], dofs_1[j], transformation(i, jj));
}
// Continue with next dof.
continue;
}
// We are left with an "identity constraint".
// Get the correct dof index on face_1 respecting the given
// orientation:
const unsigned int j = cell_to_rotated_face_index[fe.face_to_cell_index(
target, 0, face_orientation, face_flip, face_rotation)];
auto dof_left = dofs_1[j];
auto dof_right = dofs_2[i];
// If dof_left is already constrained, or dof_left < dof_right we
// flip the order to ensure that dofs are constrained in a stable
// manner on different MPI processes.
if (affine_constraints.is_constrained(dof_left) ||
(dof_left < dof_right &&
!affine_constraints.is_constrained(dof_right))) {
std::swap(dof_left, dof_right);
constraint_factor = 1. / constraint_factor;
}
// Next, we try to enter the constraint
// dof_left = constraint_factor * dof_right;
// If both degrees of freedom are constrained, there is nothing we
// can do. Simply continue with the next dof.
if (affine_constraints.is_constrained(dof_left) &&
affine_constraints.is_constrained(dof_right))
continue;
// We have to be careful that adding the current identity
// constraint does not create a constraint cycle. Thus, check for
// a dependency cycle:
bool constraints_are_cyclic = true;
number cycle_constraint_factor = constraint_factor;
for (auto test_dof = dof_right; test_dof != dof_left;) {
if (!affine_constraints.is_constrained(test_dof)) {
constraints_are_cyclic = false;
break;
}
const auto& constraint_entries =
*affine_constraints.get_constraint_entries(test_dof);
if (constraint_entries.size() == 1) {
test_dof = constraint_entries[0].first;
cycle_constraint_factor *= constraint_entries[0].second;
}
else {
constraints_are_cyclic = false;
break;
}
}
// In case of a dependency cycle we, either
// - do nothing if cycle_constraint_factor == 1. In this case all
// degrees
// of freedom are already periodically constrained,
// - otherwise, force all dofs to zero (by setting dof_left to
// zero). The reasoning behind this is the fact that
// cycle_constraint_factor != 1 occurs in situations such as
// x1 == x2 and x2 == -1. * x1. This system is only solved by
// x_1 = x_2 = 0.
if (constraints_are_cyclic) {
if (std::abs(cycle_constraint_factor - 1.) > eps)
affine_constraints.add_line(dof_left);
}
else {
affine_constraints.add_line(dof_left);
affine_constraints.add_entry(
dof_left, dof_right, constraint_factor);
// My addition for inhomogeneity
Point<dim> dof_left_support {dofs_to_supports[dof_left]};
Point<dim> dof_right_support {dofs_to_supports[dof_right]};
Point<dim> diff {dof_right_support - dof_left_support};
Vector<double> offset(3);
// This could be made more general
Point<dim> ref {
abs(diff[0]), dof_left_support[1], dof_left_support[2]};
boundary_function.vector_value(ref, offset);
affine_constraints.set_inhomogeneity(
dof_left, offset[component_index]);
// The number 1e10 in the assert below is arbitrary. If the
// absolute value of constraint_factor is too large, then probably
// the absolute value of periodicity_factor is too large or too
// small. This would be equivalent to an evanescent wave that has
// a very small wavelength. A quick calculation shows that if
// |periodicity_factor| > 1e10 -> |np.exp(ikd)|> 1e10, therefore k
// is imaginary (evanescent wave) and the evanescent wavelength is
// 0.27 times smaller than the dimension of the structure,
// lambda=((2*pi)/log(1e10))*d. Imaginary wavenumbers can be
// interesting in some cases
// (https://doi.org/10.1103/PhysRevA.94.033813).In order to
// implement the case of in which the wavevector can be imaginary
// it would be necessary to rewrite this function and the dof
// ordering method should be modified.
// Let's take the following constraint a*x1 + b*x2 = 0. You could
// just always pick x1 = b/a*x2, but in practice this is not so
// stable if a could be a small number -- intended to be zero, but
// just very small due to roundoff. Of course, constraining x2 in
// terms of x1 has the same problem. So one chooses x1 = b/a*x2 if
// |b|<|a|, and x2 = a/b*x1 if |a|<|b|.
Assert(std::abs(constraint_factor) < 1e10,
ExcMessage("The periodicity constraint is too large. The "
"parameter periodicity_factor might be too large "
"or too small."));
}
} /* for dofs_per_face */
}
template <int dim, int spacedim>
FullMatrix<double> compute_transformation(
const FiniteElement<dim, spacedim>& fe,
const FullMatrix<double>& matrix,
const std::vector<unsigned int>& first_vector_components) {
// TODO: the implementation makes the assumption that all faces have the
// same number of dofs
// AssertDimension(fe.n_unique_faces(), 1);
// const unsigned int face_no = 0;
Assert(matrix.m() == matrix.n(), ExcInternalError());
const unsigned int n_dofs_per_face = fe.n_dofs_per_face();
if (matrix.m() == n_dofs_per_face) {
// In case of m == n == n_dofs_per_face the supplied matrix is already
// an interpolation matrix, so we use it directly:
return matrix;
}
if (first_vector_components.empty() && matrix.m() == 0) {
// Just the identity matrix in case no rotation is specified:
return IdentityMatrix(n_dofs_per_face);
}
// The matrix describes a rotation and we have to build a transformation
// matrix, we assume that for a 0* rotation we would have to build the
// identity matrix
Assert(matrix.m() == spacedim, ExcInternalError())
Quadrature<dim - 1>
quadrature(fe.get_unit_face_support_points());
// have an array that stores the location of each vector-dof tuple we want
// to rotate.
using DoFTuple = std::array<unsigned int, spacedim>;
// start with a pristine interpolation matrix...
FullMatrix<double> transformation = IdentityMatrix(n_dofs_per_face);
for (unsigned int i = 0; i < n_dofs_per_face; ++i) {
std::vector<unsigned int>::const_iterator comp_it = std::find(
first_vector_components.begin(),
first_vector_components.end(),
fe.face_system_to_component_index(i).first);
if (comp_it != first_vector_components.end()) {
const unsigned int first_vector_component = *comp_it;
// find corresponding other components of vector
DoFTuple vector_dofs;
vector_dofs[0] = i;
unsigned int n_found = 1;
Assert(*comp_it + spacedim <= fe.n_components(),
ExcMessage("Error: the finite element does not have enough "
"components "
"to define rotated periodic boundaries."));
for (unsigned int k = 0; k < n_dofs_per_face; ++k)
if ((k != i) && (quadrature.point(k) == quadrature.point(i)) &&
(fe.face_system_to_component_index(k).first >=
first_vector_component) &&
(fe.face_system_to_component_index(k).first <
first_vector_component + spacedim)) {
vector_dofs
[fe.face_system_to_component_index(k).first -
first_vector_component] = k;
n_found++;
if (n_found == dim)
break;
}
// ... and rotate all dofs belonging to vector valued components
// that are selected by first_vector_components:
for (int i = 0; i < spacedim; ++i) {
transformation[vector_dofs[i]][vector_dofs[i]] = 0.;
for (int j = 0; j < spacedim; ++j)
transformation[vector_dofs[i]][vector_dofs[j]] =
matrix[i][j];
}
}
}
return transformation;
}
template <typename FaceIterator, typename number, int dim>
void make_periodicity_constraints(
const FaceIterator& face_1,
const typename identity<FaceIterator>::type& face_2,
AffineConstraints<number>& constraints,
std::vector<Point<dim>> dofs_to_supports,
Function<dim>& boundary_function,
const ComponentMask& component_mask = ComponentMask(),
const bool face_orientation = true,
const bool face_flip = false,
const bool face_rotation = false,
const FullMatrix<double>& matrix = FullMatrix<double>(),
const std::vector<unsigned int>& first_vector_components =
std::vector<unsigned int>(),
const number periodicity_factor = 1.);
template <typename FaceIterator, typename number, int dim>
void make_periodicity_constraints(
const FaceIterator& face_1,
const typename identity<FaceIterator>::type& face_2,
AffineConstraints<number>& affine_constraints,
std::vector<Point<dim>> dofs_to_supports,
Function<dim>& boundary_function,
const ComponentMask& component_mask,
const bool face_orientation,
const bool face_flip,
const bool face_rotation,
const FullMatrix<double>& matrix,
const std::vector<unsigned int>& first_vector_components,
const number periodicity_factor) {
// TODO: the implementation makes the assumption that all faces have the
// same number of dofs
// AssertDimension(
// face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(),
// 1);
// AssertDimension(
// face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(),
// 1);
// const unsigned int face_no = 0;
static const int spacedim = FaceIterator::AccessorType::space_dimension;
Assert((dim != 1) || (face_orientation == true && face_flip == false &&
face_rotation == false),
ExcMessage("The supplied orientation "
"(face_orientation, face_flip, face_rotation) "
"is invalid for 1D"));
Assert((dim != 2) || (face_orientation == true && face_rotation == false),
ExcMessage("The supplied orientation "
"(face_orientation, face_flip, face_rotation) "
"is invalid for 2D"));
Assert(face_1 != face_2,
ExcMessage("face_1 and face_2 are equal! Cannot constrain DoFs "
"on the very same face"));
Assert(face_1->at_boundary() && face_2->at_boundary(),
ExcMessage("Faces for periodicity constraints must be on the "
"boundary"));
Assert(matrix.m() == matrix.n(),
ExcMessage("The supplied (rotation or interpolation) matrix must "
"be a square matrix"));
Assert(first_vector_components.empty() || matrix.m() == spacedim,
ExcMessage("first_vector_components is nonempty, so matrix must "
"be a rotation matrix exactly of size spacedim"));
#ifdef DEBUG
if (!face_1->has_children()) {
Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
const unsigned int n_dofs_per_face =
face_1->get_fe(face_1->nth_active_fe_index(0))
.n_dofs_per_face();
Assert(matrix.m() == 0 ||
(first_vector_components.empty() &&
matrix.m() == n_dofs_per_face) ||
(!first_vector_components.empty() &&
matrix.m() == spacedim),
ExcMessage(
"The matrix must have either size 0 or spacedim "
"(if first_vector_components is nonempty) "
"or the size must be equal to the # of DoFs on the face "
"(if first_vector_components is empty)."));
}
if (!face_2->has_children()) {
Assert(face_2->n_active_fe_indices() == 1, ExcInternalError());
const unsigned int n_dofs_per_face =
face_2->get_fe(face_2->nth_active_fe_index(0))
.n_dofs_per_face();
Assert(matrix.m() == 0 ||
(first_vector_components.empty() &&
matrix.m() == n_dofs_per_face) ||
(!first_vector_components.empty() &&
matrix.m() == spacedim),
ExcMessage(
"The matrix must have either size 0 or spacedim "
"(if first_vector_components is nonempty) "
"or the size must be equal to the # of DoFs on the face "
"(if first_vector_components is empty)."));
}
#endif
// A lookup table on how to go through the child faces depending on the
// orientation:
static const int lookup_table_2d[2][2] = {
// flip:
{0, 1}, // false
{1, 0}, // true
};
static const int lookup_table_3d[2][2][2][4] = {
// orientation flip rotation
{
{
{0, 2, 1, 3}, // false false false
{2, 3, 0, 1}, // false false true
},
{
{3, 1, 2, 0}, // false true false
{1, 0, 3, 2}, // false true true
},
},
{
{
{0, 1, 2, 3}, // true false false
{1, 3, 0, 2}, // true false true
},
{
{3, 2, 1, 0}, // true true false
{2, 0, 3, 1}, // true true true
},
},
};
if (face_1->has_children() && face_2->has_children()) {
// In the case that both faces have children, we loop over all children
// and apply make_periodicty_constrains recursively:
Assert(face_1->n_children() ==
GeometryInfo<dim>::max_children_per_face &&
face_2->n_children() ==
GeometryInfo<dim>::max_children_per_face,
ExcNotImplemented());
for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face;
++i) {
// Lookup the index for the second face
unsigned int j;
switch (dim) {
case 2:
j = lookup_table_2d[face_flip][i];
break;
case 3:
j = lookup_table_3d[face_orientation][face_flip][face_rotation]
[i];
break;
default:
AssertThrow(false, ExcNotImplemented());
}
make_periodicity_constraints(
face_1->child(i),
face_2->child(j),
affine_constraints,
dofs_to_supports,
boundary_function,
component_mask,
face_orientation,
face_flip,
face_rotation,
matrix,
first_vector_components,
periodicity_factor);
}
}
else {
// Otherwise at least one of the two faces is active and we need to do
// some work and enter the constraints!
// The finite element that matters is the one on the active face:
const FiniteElement<dim, spacedim>& fe =
face_1->has_children() ?
face_2->get_fe(face_2->nth_active_fe_index(0)) :
face_1->get_fe(face_1->nth_active_fe_index(0));
const unsigned int n_dofs_per_face = fe.n_dofs_per_face();
// Sometimes we just have nothing to do (for all finite elements, or
// systems which accidentally don't have any dofs on the boundary).
if (n_dofs_per_face == 0)
return;
const FullMatrix<double> transformation =
compute_transformation(fe, matrix, first_vector_components);
if (!face_2->has_children()) {
// Performance hack: We do not need to compute an inverse if the
// matrix is the identity matrix.
if (first_vector_components.empty() && matrix.m() == 0) {
set_periodicity_constraints(
face_2,
face_1,
transformation,
affine_constraints,
dofs_to_supports,
boundary_function,
component_mask,
face_orientation,
face_flip,
face_rotation,
periodicity_factor);
}
else {
FullMatrix<double> inverse(transformation.m());
inverse.invert(transformation);
set_periodicity_constraints(
face_2,
face_1,
inverse,
affine_constraints,
dofs_to_supports,
boundary_function,
component_mask,
face_orientation,
face_flip,
face_rotation,
periodicity_factor);
}
}
else {
Assert(!face_1->has_children(), ExcInternalError());
// Important note:
// In 3D we have to take care of the fact that face_rotation gives
// the relative rotation of face_1 to face_2, i.e. we have to invert
// the rotation when constraining face_2 to face_1. Therefore
// face_flip has to be toggled if face_rotation is true: In case of
// inverted orientation, nothing has to be done.
set_periodicity_constraints(
face_1,
face_2,
transformation,
affine_constraints,
dofs_to_supports,
boundary_function,
component_mask,
face_orientation,
face_orientation ? face_rotation ^ face_flip : face_flip,
face_rotation,
periodicity_factor);
}
}
}
template <int dim, int spacedim, typename number>
void make_periodicity_constraints(
const std::vector<GridTools::PeriodicFacePair<
typename DoFHandler<dim, spacedim>::cell_iterator>>&
periodic_faces,
AffineConstraints<number>& constraints,
std::vector<Point<dim>> dofs_to_supports,
Function<dim>& boundary_function,
const ComponentMask& component_mask = ComponentMask(),
const std::vector<unsigned int>& first_vector_components =
std::vector<unsigned int>(),
const number periodicity_factor = 1.);
template <int dim, int spacedim, typename number>
void make_periodicity_constraints(
const std::vector<GridTools::PeriodicFacePair<
typename DoFHandler<dim, spacedim>::cell_iterator>>&
periodic_faces,
AffineConstraints<number>& constraints,
std::vector<Point<dim>> dofs_to_supports,
Function<dim>& boundary_function,
const ComponentMask& component_mask,
const std::vector<unsigned int>& first_vector_components,
const number periodicity_factor) {
// Loop over all periodic faces...
for (auto& pair: periodic_faces) {
using FaceIterator = typename DoFHandler<dim, spacedim>::face_iterator;
const FaceIterator face_1 = pair.cell[0]->face(pair.face_idx[0]);
const FaceIterator face_2 = pair.cell[1]->face(pair.face_idx[1]);
Assert(face_1->at_boundary() && face_2->at_boundary(),
ExcInternalError());
Assert(face_1 != face_2, ExcInternalError());
// ... and apply the low level make_periodicity_constraints function to
// every matching pair:
make_periodicity_constraints(
face_1,
face_2,
constraints,
dofs_to_supports,
boundary_function,
component_mask,
pair.orientation[0],
pair.orientation[1],
pair.orientation[2],
pair.matrix,
first_vector_components,
periodicity_factor);
}
}
#endif
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/parsed_function.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/utilities.h>
#include <deal.II/differentiation/ad.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/data_out_faces.h>
#include <deal.II/numerics/data_postprocessor.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/vector_tools.h>
#include "modded_periodic_functions.h"
using namespace dealii;
template <int dim>
struct Params {
unsigned int global_refinements;
unsigned int adaptive_refinements;
double length;
double width;
unsigned int num_boundary_stages;
unsigned int starting_stage;
std::vector<unsigned int> num_gamma_iters;
// There is probably a better way to do this
static const unsigned int max_stages {5};
std::vector<std::shared_ptr<Functions::ParsedFunction<dim>>>
boundary_functions;
double E;
double nu;
unsigned int max_n_line_searches;
double alpha_factor;
double min_alpha;
double alpha_check_factor;
double nonlinear_tol;
std::string linear_solver;
unsigned int max_linear_iters;
double linear_tol;
double precon_relaxation_param;
bool load_from_checkpoint;
std::string input_prefix;
std::string input_checkpoint;
std::string output_prefix;
unsigned int gamma_precision;
Params();
static void declare_parameters(ParameterHandler& prm);
void parse_parameters(ParameterHandler& prm);
};
template <int dim>
Params<dim>::Params() {
for (unsigned int i {0}; i != max_stages; i++) {
boundary_functions.push_back(
std::make_shared<Functions::ParsedFunction<dim>>(dim));
}
ParameterHandler prm;
declare_parameters(prm);
prm.parse_input(std::cin);
parse_parameters(prm);
}
template <int dim>
void Params<dim>::declare_parameters(ParameterHandler& prm) {
prm.enter_subsection("Mesh and geometry");
{
prm.declare_entry(
"Number of initial refinements",
"2",
Patterns::Integer(0),
"Number of initial mesh refinements");
prm.declare_entry(
"Number of final refinements",
"0",
Patterns::Integer(0),
"Number of final mesh refinements");
prm.declare_entry(
"Beam length", "20", Patterns::Double(0), "Length of beam");
prm.declare_entry(
"Beam width", "2", Patterns::Double(0), "Width of beam");
}
prm.leave_subsection();
prm.enter_subsection("Boundary conditions");
{
prm.declare_entry(
"Number of boundary stages",
"1",
Patterns::Integer(1),
"Number of boundary stages");
prm.declare_entry(
"Starting stage",
"0",
Patterns::Integer(0),
"Boundary stage to start calculations on");
}
prm.leave_subsection();
for (unsigned int i {0}; i != max_stages + 1; i++) {
prm.enter_subsection("Boundary function stage " + std::to_string(i));
{
prm.declare_entry(
"Number of boundary increments",
"1",
Patterns::Integer(1),
"Number of boundary increments");
Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
}
prm.leave_subsection();
}
prm.enter_subsection("Physical constants");
{
prm.declare_entry(
"Young's modulus",
"2.2e7",
Patterns::Double(0),
"Young's modulus");
prm.declare_entry(
"Poisson's ratio",
"0.3",
Patterns::Double(0),
"Poisson's ratio");
}
prm.leave_subsection();
prm.enter_subsection("Solver");
{
prm.declare_entry(
"Maximum number of line searches",
"10",
Patterns::Integer(0),
"Maximum number of line searches");
prm.declare_entry(
"Minimum alpha",
"1e-5",
Patterns::Double(0),
"Smallest alpha to try before moving on");
prm.declare_entry(
"Alpha factor",
"0.5",
Patterns::Double(0),
"Factor to multiply alpha by when trying smaller step");
prm.declare_entry(
"Alpha check factor",
"0.5",
Patterns::Double(0),
"Larger values require bigger residual decrease for acceptable "
"alpha");
prm.declare_entry(
"Nonlinear tolerance",
"1e-12",
Patterns::Double(0),
"Tolerance for nonlinear solver");
prm.declare_entry(
"Linear solver",
"direct",
Patterns::Anything(),
"Type of linear solver (direct or iterative)");
prm.declare_entry(
"Maximum number of linear solver iterations",
"10000",
Patterns::Integer(1),
"Maximum number of linear solver iterations");
prm.declare_entry(
"Linear tolerance",
"1e-12",
Patterns::Double(0),
"Tolerance for linear solver");
prm.declare_entry(
"Preconditioner relaxation parameter",
"1.0",
Patterns::Double(0),
"Preconditioner relaxation parameter");
}
prm.leave_subsection();
prm.enter_subsection("Input and output");
{
prm.declare_entry(
"Load from checkpoint",
"false",
Patterns::Bool(),
"Begin calculation from a checkpoint");
prm.declare_entry(
"Input filename prefix",
"",
Patterns::Anything(),
"Prefix of the output filename");
prm.declare_entry(
"Input checkpoint",
"",
Patterns::Anything(),
"Checkpoint of the input file");
prm.declare_entry(
"Output filename prefix",
"",
Patterns::Anything(),
"Prefix of the output filename");
prm.declare_entry(
"Gamma precision",
"1",
Patterns::Integer(0),
"Number of digits to include in gamma string");
}
prm.leave_subsection();
}
template <int dim>
void Params<dim>::parse_parameters(ParameterHandler& prm) {
prm.enter_subsection("Mesh and geometry");
{
global_refinements = prm.get_integer("Number of initial refinements");
adaptive_refinements = prm.get_integer("Number of final refinements");
length = prm.get_double("Beam length");
width = prm.get_double("Beam width");
}
prm.leave_subsection();
prm.enter_subsection("Boundary conditions");
{
num_boundary_stages = prm.get_integer("Number of boundary stages");
starting_stage = prm.get_integer("Starting stage");
}
prm.leave_subsection();
for (unsigned int i {0}; i != num_boundary_stages + 1; i++) {
prm.enter_subsection("Boundary function stage " + std::to_string(i));
{
num_gamma_iters.push_back(
prm.get_integer("Number of boundary increments"));
boundary_functions[i]->parse_parameters(prm);
}
prm.leave_subsection();
}
prm.enter_subsection("Physical constants");
{
E = prm.get_double("Young's modulus");
nu = prm.get_double("Poisson's ratio");
}
prm.leave_subsection();
prm.enter_subsection("Solver");
{
max_n_line_searches =
prm.get_integer("Maximum number of line searches");
alpha_factor = prm.get_double("Alpha factor");
alpha_check_factor = prm.get_double("Alpha check factor");
min_alpha = prm.get_double("Minimum alpha");
nonlinear_tol = prm.get_double("Nonlinear tolerance");
linear_solver = prm.get("Linear solver");
max_linear_iters =
prm.get_integer("Maximum number of linear solver iterations");
linear_tol = prm.get_double("Linear tolerance");
precon_relaxation_param =
prm.get_double("Preconditioner relaxation parameter");
}
prm.leave_subsection();
prm.enter_subsection("Input and output");
{
load_from_checkpoint = prm.get_bool("Load from checkpoint");
input_prefix = prm.get("Input filename prefix");
input_checkpoint = prm.get("Input checkpoint");
output_prefix = prm.get("Output filename prefix");
gamma_precision = prm.get_integer("Gamma precision");
}
prm.leave_subsection();
}
template <int dim>
class ComposedFunction: public Function<dim> {
public:
ComposedFunction();
void vector_value(const Point<dim>& p, Vector<double>& values)
const override;
std::shared_ptr<Function<dim>> f1;
std::shared_ptr<Function<dim>> f2;
double gamma {1};
};
template <int dim>
ComposedFunction<dim>::ComposedFunction() {
f1 = std::make_shared<Functions::ZeroFunction<dim>>(3);
f2 = std::make_shared<Functions::ZeroFunction<dim>>(3);
}
template <int dim>
void ComposedFunction<dim>::vector_value(
const Point<dim>& p,
Vector<double>& values) const {
Vector<double> values_1(3);
f1->vector_value(p, values_1);
Vector<double> values_2(3);
f2->vector_value(p, values_2);
values.add(1 - gamma, values_1, gamma, values_2);
}
template <int dim>
class SolveRing {
public:
SolveRing(Params<dim>& prms);
void run();
private:
void make_grid();
void initiate_system();
void setup_constraints(unsigned int stage_i);
void update_constraints();
void setup_sparsity_pattern();
void assemble(const bool initial_step, const bool assemble_matrix);
void assemble_system(const bool initial_step);
void assemble_rhs(const bool initial_step);
void solve(const bool initial_step);
double calc_residual_norm();
void newton_iteration(bool first_step, const std::string checkpoint);
void output_grid() const;
void output_checkpoint(const std::string checkpoint) const;
void output_results(const std::string checkpoint) const;
void load_checkpoint(const std::string checkpoint);
std::string format_gamma();
void update_boundary_function_gamma();
void update_boundary_function_stage(unsigned int stage_i);
void refine_mesh(unsigned int stage_i);
void center_solution_on_mean();
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> zero_constraints;
AffineConstraints<double> nonzero_constraints;
Functions::ZeroFunction<dim> zero_function;
ComposedFunction<dim> boundary_function;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> present_solution;
Vector<double> newton_update;
Vector<double> system_rhs;
Vector<double> evaluation_point;
double present_energy;
std::vector<Point<dim>> dofs_to_supports;
std::vector<GridTools::PeriodicFacePair<
typename DoFHandler<dim, dim>::cell_iterator>>
face_pairs;
double gamma;
Params<dim>& prms;
double lambda;
double mu;
};
template <int dim>
SolveRing<dim>::SolveRing(Params<dim>& prms):
fe(FE_Q<dim>(1), dim),
dof_handler(triangulation),
zero_function {3},
prms {prms} {
mu = prms.E / (2 * (1 + prms.nu));
lambda = prms.E * prms.nu / ((1 + prms.nu) * (1 - 2 * prms.nu));
}
template <int dim>
void SolveRing<dim>::make_grid() {
const Point<dim>& origin {0, 0, 0};
const Point<dim>& size {prms.length, prms.width, prms.width};
GridGenerator::hyper_rectangle(triangulation, origin, size);
for (auto& face: triangulation.active_face_iterators()) {
if (std::fabs(face->center()(1) - prms.width / 2) < 1e-12) {
if (std::fabs(face->center()(0)) < 1e-12) {
face->set_boundary_id(1);
}
else if (std::fabs(face->center()(0) - prms.length) < 1e-12) {
face->set_boundary_id(2);
}
}
}
triangulation.refine_global(prms.global_refinements);
}
template <int dim>
void SolveRing<dim>::initiate_system() {
dof_handler.distribute_dofs(fe);
present_solution.reinit(dof_handler.n_dofs());
newton_update.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void SolveRing<dim>::setup_constraints(unsigned int stage_i) {
const MappingQ1<dim> mapping;
dofs_to_supports.resize(dof_handler.n_dofs());
DoFTools::map_dofs_to_support_points<dim, dim>(
mapping, dof_handler, dofs_to_supports);
collect_periodic_faces(dof_handler, 1, 2, 0, face_pairs);
update_boundary_function_stage(stage_i);
nonzero_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints);
make_periodicity_constraints<dim, dim, double>(
face_pairs,
nonzero_constraints,
dofs_to_supports,
boundary_function);
nonzero_constraints.close();
zero_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints);
make_periodicity_constraints<dim, dim, double>(
face_pairs, zero_constraints, dofs_to_supports, zero_function);
zero_constraints.close();
cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
}
template <int dim>
void SolveRing<dim>::update_constraints() {
update_boundary_function_gamma();
nonzero_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints);
make_periodicity_constraints<dim, dim, double>(
face_pairs,
nonzero_constraints,
dofs_to_supports,
boundary_function);
nonzero_constraints.close();
}
template <int dim>
void SolveRing<dim>::setup_sparsity_pattern() {
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp, nonzero_constraints);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
}
template <int dim>
void SolveRing<dim>::assemble(
const bool initial_step,
const bool assemble_matrix) {
if (assemble_matrix) {
system_matrix = 0;
}
system_rhs = 0;
present_energy = 0;
QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(
fe,
quadrature_formula,
update_values | update_gradients | update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
const FEValuesExtractors::Vector displacements(0);
using ADHelper = Differentiation::AD::EnergyFunctional<
Differentiation::AD::NumberTypes::sacado_dfad_dfad,
double>;
using ADNumberType = typename ADHelper::ad_type;
for (const auto& cell: dof_handler.active_cell_iterators()) {
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
const unsigned int n_independent_variables = local_dof_indices.size();
ADHelper ad_helper(n_independent_variables);
ad_helper.register_dof_values(evaluation_point, local_dof_indices);
const std::vector<ADNumberType>& dof_values_ad =
ad_helper.get_sensitive_dof_values();
ADNumberType energy_ad = ADNumberType(0.0);
std::vector<Tensor<2, dim, ADNumberType>> old_solution_gradients(
fe_values.n_quadrature_points);
fe_values[displacements].get_function_gradients_from_local_dof_values(
dof_values_ad, old_solution_gradients);
for (const unsigned int q_index: fe_values.quadrature_point_indices()) {
const Tensor<2, dim, ADNumberType> grad_u {
old_solution_gradients[q_index]};
const Tensor<2, dim, ADNumberType> grad_u_T {transpose(grad_u)};
const Tensor<2, dim, ADNumberType> green_lagrange_strain_tensor {
0.5 * (grad_u + grad_u_T + grad_u * grad_u_T)};
ADNumberType t1 = lambda / 2 *
std::pow(trace(green_lagrange_strain_tensor), 2);
ADNumberType t2 = mu * trace(green_lagrange_strain_tensor *
green_lagrange_strain_tensor);
ADNumberType pi {t1 + t2};
energy_ad += pi * fe_values.JxW(q_index);
}
ad_helper.register_energy_functional(energy_ad);
present_energy += ad_helper.compute_energy();
ad_helper.compute_residual(cell_rhs);
cell_rhs *= -1.0; // RHS = - residual
if (assemble_matrix) {
ad_helper.compute_linearization(cell_matrix);
}
const AffineConstraints<double>& constraints_used =
initial_step ? nonzero_constraints : zero_constraints;
if (assemble_matrix) {
constraints_used.distribute_local_to_global(
cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs,
true);
}
else {
constraints_used.distribute_local_to_global(
cell_rhs, local_dof_indices, system_rhs);
}
}
}
template <int dim>
void SolveRing<dim>::assemble_system(const bool initial_step) {
assemble(initial_step, true);
}
template <int dim>
void SolveRing<dim>::assemble_rhs(const bool initial_step) {
assemble(initial_step, false);
}
template <int dim>
void SolveRing<dim>::solve(const bool initial_step) {
newton_update = 0;
const AffineConstraints<double>& constraints_used =
initial_step ? nonzero_constraints : zero_constraints;
if (prms.linear_solver == "iterative") {
SolverControl solver_control(prms.max_linear_iters, prms.linear_tol);
SolverGMRES<Vector<double>> solver(solver_control);
PreconditionSSOR<SparseMatrix<double>> preconditioner;
preconditioner.initialize(system_matrix, prms.precon_relaxation_param);
solver.solve(system_matrix, newton_update, system_rhs, preconditioner);
}
else if (prms.linear_solver == "direct") {
SparseDirectUMFPACK A_direct;
A_direct.initialize(system_matrix);
A_direct.vmult(newton_update, system_rhs);
}
constraints_used.distribute(newton_update);
}
template <int dim>
double SolveRing<dim>::calc_residual_norm() {
Vector<double> residual(dof_handler.n_dofs());
for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) {
if (!nonzero_constraints.is_constrained(i)) {
residual(i) = system_rhs(i);
}
}
return residual.l2_norm();
}
template <int dim>
void SolveRing<dim>::newton_iteration(
bool first_step,
const std::string checkpoint) {
bool boundary_updated {!first_step};
unsigned int line_search_n = 0;
double last_res = 1.0;
double current_res = 1.0;
while ((first_step or (current_res > prms.nonlinear_tol)) and
line_search_n < prms.max_n_line_searches) {
if (first_step) {
evaluation_point = present_solution;
assemble_system(first_step);
solve(first_step);
present_solution = newton_update;
first_step = false;
evaluation_point = present_solution;
assemble_rhs(first_step);
current_res = calc_residual_norm();
std::cout << "The residual of initial guess is " << current_res
<< std::endl;
std::cout << "The energy of initial guess is " << present_energy
<< std::endl;
last_res = current_res;
}
else {
nonzero_constraints.distribute(present_solution);
evaluation_point = present_solution;
assemble_system(first_step);
solve(first_step);
for (double alpha {1.0}; alpha > prms.min_alpha;
alpha *= prms.alpha_factor) {
evaluation_point = present_solution;
evaluation_point.add(alpha, newton_update);
assemble_rhs(first_step);
current_res = calc_residual_norm();
std::cout << " alpha: " << std::setw(10) << alpha
<< std::setw(0) << " residual: " << current_res
<< std::endl;
if (boundary_updated or
((last_res - current_res) >=
(alpha * last_res * prms.alpha_check_factor))) {
boundary_updated = false;
break;
}
}
present_solution = evaluation_point;
std::cout << " number of line searches: " << line_search_n
<< " residual: " << current_res
<< " energy: " << present_energy << std::endl;
last_res = current_res;
++line_search_n;
center_solution_on_mean();
}
}
output_results(checkpoint);
}
template <int dim>
void SolveRing<dim>::center_solution_on_mean() {
Vector<double> mean_u(3);
// This would need to be changed to only count each global dof once
/*QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe, quadrature_formula, update_default);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
for (const auto& cell: dof_handler.active_cell_iterators()) {
fe_values.reinit(cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
for (const unsigned int i: fe_values.dof_indices()) {
const unsigned int component_i {
fe.system_to_component_index(i).first};
const unsigned int global_dof_i {local_dof_indices[i]};
mean_u[component_i] += present_solution[global_dof_i];
}
}
mean_u /= static_cast<double>(dof_handler.n_dofs()) / 3;
for (const auto& cell: dof_handler.active_cell_iterators()) {
fe_values.reinit(cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
for (const unsigned int i: fe_values.dof_indices()) {
const unsigned int component_i {
fe.system_to_component_index(i).first};
const unsigned int global_dof_i {local_dof_indices[i]};
present_solution[global_dof_i] -= mean_u[component_i];
}
}*/
for (unsigned int i {0}; i != dof_handler.n_dofs(); i++) {
const unsigned int component_i {i % dim};
mean_u[component_i] += present_solution[i];
}
mean_u /= static_cast<double>(dof_handler.n_dofs()) / 3;
for (unsigned int i {0}; i != dof_handler.n_dofs(); i++) {
const unsigned int component_i {i % dim};
present_solution[i] -= mean_u[component_i];
}
}
template <int dim>
void SolveRing<dim>::output_grid() const {
GridOut grid_out;
std::ofstream grid_output {prms.output_prefix + "_mesh.vtk"};
grid_out.write_vtk(triangulation, grid_output);
}
template <int dim>
void SolveRing<dim>::output_checkpoint(const std::string checkpoint) const {
std::ofstream solution_out {
prms.output_prefix + "_displacement_" + checkpoint + ".ar"};
boost::archive::text_oarchive solution_ar {solution_out};
present_solution.save(solution_ar, 0);
std::ofstream triangulation_out {
prms.output_prefix + "_triangulation_" + checkpoint + ".ar"};
boost::archive::text_oarchive triangulation_ar {triangulation_out};
triangulation.save(triangulation_ar, 0);
std::ofstream dof_handler_out {
prms.output_prefix + "_dof_handler_" + checkpoint + ".ar"};
boost::archive::text_oarchive dof_handler_ar {dof_handler_out};
dof_handler.save(dof_handler_ar, 0);
}
template <int dim>
void SolveRing<dim>::output_results(const std::string checkpoint) const {
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(
dim,
DataComponentInterpretation::component_is_part_of_vector);
std::vector<std::string> solution_names(dim, "displacement");
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(
present_solution,
solution_names,
DataOut<dim>::type_dof_data,
data_component_interpretation);
data_out.build_patches();
std::ofstream data_output(
prms.output_prefix + "_cells_" + checkpoint + ".vtk");
data_out.write_vtk(data_output);
}
template <int dim>
void SolveRing<dim>::load_checkpoint(const std::string checkpoint) {
std::ifstream solution_inp {
prms.input_prefix + "_displacement_" + checkpoint + ".ar"};
boost::archive::text_iarchive solution_ar {solution_inp};
present_solution.load(solution_ar, 0);
std::ifstream triangulation_inp {
prms.input_prefix + "_triangulation_" + checkpoint + ".ar"};
boost::archive::text_iarchive triangulation_ar {triangulation_inp};
triangulation.load(triangulation_ar, 0);
dof_handler.distribute_dofs(fe);
std::ifstream dof_handler_inp {
prms.input_prefix + "_dof_handler_" + checkpoint + ".ar"};
boost::archive::text_iarchive dof_handler_ar {dof_handler_inp};
dof_handler.load(dof_handler_ar, 0);
newton_update.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
std::string SolveRing<dim>::format_gamma() {
std::ostringstream stream_obj;
stream_obj << std::fixed;
stream_obj << std::setprecision(prms.gamma_precision);
stream_obj << gamma;
return stream_obj.str();
}
template <int dim>
void SolveRing<dim>::update_boundary_function_gamma() {
boundary_function.gamma = gamma;
}
template <int dim>
void SolveRing<dim>::update_boundary_function_stage(unsigned int stage_i) {
boundary_function.f1 = prms.boundary_functions[stage_i];
boundary_function.f2 = prms.boundary_functions[stage_i + 1];
}
template <int dim>
void SolveRing<dim>::refine_mesh(unsigned int stage_i) {
/*Vector<float>
estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(
dof_handler,
QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim>*>(),
present_solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_number(
triangulation, estimated_error_per_cell, 0.3, 0.0);
triangulation.prepare_coarsening_and_refinement();
SolutionTransfer<dim, Vector<double>> soltrans(dof_handler);
soltrans.prepare_for_coarsening_and_refinement(present_solution);
triangulation.execute_coarsening_and_refinement();*/
Vector<double> cells_to_refine(triangulation.n_active_cells());
cells_to_refine = 1;
GridRefinement::refine(triangulation, cells_to_refine, 0);
triangulation.prepare_coarsening_and_refinement();
SolutionTransfer<dim, Vector<double>> soltrans(dof_handler);
soltrans.prepare_for_pure_refinement();
triangulation.execute_coarsening_and_refinement();
dof_handler.distribute_dofs(fe);
Vector<double> interpolated_solution(dof_handler.n_dofs());
soltrans.refine_interpolate(present_solution, interpolated_solution);
present_solution.reinit(dof_handler.n_dofs());
newton_update.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
setup_constraints(stage_i);
setup_sparsity_pattern();
present_solution = interpolated_solution;
}
template <int dim>
void SolveRing<dim>::run() {
make_grid();
output_grid();
std::string checkpoint;
bool first_step {true};
if (prms.load_from_checkpoint) {
first_step = false;
load_checkpoint(prms.input_checkpoint);
}
else {
initiate_system();
}
setup_constraints(prms.starting_stage);
setup_sparsity_pattern();
for (unsigned int stage_i {prms.starting_stage};
stage_i != prms.num_boundary_stages;
stage_i++) {
update_boundary_function_stage(stage_i);
unsigned int num_gamma_iters {prms.num_gamma_iters[stage_i]};
for (unsigned int i {0}; i != num_gamma_iters; i++) {
gamma = static_cast<double>((i + 1)) / num_gamma_iters;
std::string gamma_formatted {format_gamma()};
checkpoint = std::to_string(stage_i + 1) + "-" + gamma_formatted;
cout << "Stage " << std::to_string(stage_i + 1) << ", gamma "
<< gamma_formatted << std::endl;
update_constraints();
newton_iteration(first_step, checkpoint);
output_checkpoint(checkpoint);
first_step = false;
}
}
checkpoint = "refine-0";
output_checkpoint(checkpoint);
output_results(checkpoint);
for (unsigned int i {0}; i != prms.adaptive_refinements; i++) {
cout << "Grid refinement " << std::to_string(i + 1) << std::endl;
checkpoint = "refine-" + std::to_string(i + 1);
refine_mesh(prms.num_boundary_stages - 1);
newton_iteration(first_step, checkpoint);
}
checkpoint = "final";
output_checkpoint(checkpoint);
output_results(checkpoint);
checkpoint = "moved-mesh";
}
int main() {
deallog.depth_console(0);
Params<3> prms {};
SolveRing<3> ring_solver {prms};
ring_solver.run();
return 0;
}
