Dear all,
*Background :* I am currently trying to implement a diffuse domain
formulation to solve the poisson equation whilst penalizing the jumps in
gradient in very similar way to what step-85 showcases. The nuance relative
to the mentioned tutorial stems the need for AMR close to the boundary.
Having worked through this problem without neglecting the outside of the
domain (contrary to what is shown in step-85 with the use of FENothing in
the FECollection object for all the elements for which the level set is
positive on all vertices), I found myself wanting to implement this feature
again.
*Issue :* Taking good care to assigne FENothing only to the cells that
extend beyond an arbitrary threshold to avoid computing the face values on
a cell who's neighbor is neglected, I got the following error when
computing the flux sparsity pattern. Note that hanging nodes are taken into
account. Note that this issue arise only if cells are flagged for ghost
penalty and if elements (outside the domain of interest) are attributed the
FENothing finite element.
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
const unsigned int n_components =
fe_collection.n_components();
Table<2, DoFTools::Coupling> cell_coupling(n_components, n_components);
Table<2, DoFTools::Coupling> face_coupling(n_components, n_components);
cell_coupling[0][0] = DoFTools::always;
face_coupling[0][0] = DoFTools::always;
const bool keep_constrained_dofs = true;
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
constraints.close();
DoFTools::make_flux_sparsity_pattern(dof_handler,
dsp,
constraints,
keep_constrained_dofs,
cell_coupling,
face_coupling,
numbers::invalid_subdomain_id);
-------------------------------------------------------------
The violated condition was:
matrix_values->column() == column
Additional information:
You are trying to access the matrix entry with index <14,142>, but
this entry does not exist in the sparsity pattern of this matrix.
-------------------------------------------------------------
Having tried the following sister make_flux_sparsity_pattern() function
call without explicitly passing the cell and face couplings, I get this
error:
-------------------------------------------------------------
An error occurred in line <4467> of file
</home/mh2294/codes/dealii-9.5.2/include/deal.II/lac/affine_constraints.templates.h>
in function
void
dealii::AffineConstraints<number>::add_entries_local_to_global(const
std::vector<unsigned int>&, const dealii::AffineConstraints<number>&, const
std::vector<unsigned int>&, dealii::SparsityPatternBase&, bool, const
dealii::Table<2, bool>&) const [with number = double]
The violated condition was:
false
Additional information:
You are trying to use functionality in deal.II that is currently not
implemented. In many cases, this indicates that there simply didn't
appear much of a need for it, or that the author of the original code
did not have the time to implement a particular case. If you hit this
exception, it is therefore worth the time to look into the code to
find out whether you may be able to implement the missing
functionality. If you do, please consider providing a patch to the
deal.II development sources (see the deal.II website on how to
contribute).
------------------------------------------------------------
Implying that the both sparsity patterns aren't the same although
specifying should result in a full coupling in the case of scalar valued
problems.
cell_coupling[0][0] = DoFTools::always;
face_coupling[0][0] = DoFTools::always;
*Questions :*
1. Is there a way to deal make a correct sparsity pattern accounting for
hanging nodes, FENothing elements and coupling of faces that have ghost
penalty or is this really a feature that isn't implemented in deal.ii ?
2. Should I assemble the matrix without taking care of the constraints and
apply them afterwards ?
3. Why do the sparsity patterns differ from one call to another ?
I have attached the slightly modified step-85 code I am talking about.
I hope this question isn't redundant,
Thank you for your help,
Matthias
--
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/fb5f6087-b96c-47f1-b1d2-36444f3192b3n%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2021 - 2023 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*/
// @sect3{Include files}
// The first include files have all been treated in previous examples.
#include <deal.II/base/function.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_interface_values.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/q_collection.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/constrained_linear_operator.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_gmres.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <fstream>
#include <vector>
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <tuple>
#include <string>
#include <deal.II/base/function_signed_distance.h>
#include <deal.II/non_matching/fe_immersed_values.h>
#include <deal.II/non_matching/fe_values.h>
#include <deal.II/non_matching/mesh_classifier.h>
namespace Step85
{
using namespace dealii;
template <int dim, typename SignedDistanceType>
class PhaseFieldFunction : public Function<dim> {
public:
PhaseFieldFunction(double eps, const SignedDistanceType &signed_distance)
: Function<dim>(), epsilon(eps), signed_distance(signed_distance) {}
void set_epsilon(double eps){
epsilon = eps;
}
virtual double value(const Point<dim> &p, const unsigned int component = 0) const override {
double signed_distance_value = signed_distance.value(p);
return .5 * (1 - std::tanh(3. * signed_distance_value / epsilon));
}
Tensor<1, dim> gradient(const Point<dim> &p, const unsigned int component = 0) const override {
double signed_distance_value = signed_distance.value(p);
Tensor<1, dim> grad_signed_distance = signed_distance.gradient(p);
Tensor<1, dim> grad;
double factor = -1.5 / epsilon * (1 / std::cosh(3 * signed_distance_value / epsilon)) * (1 / std::cosh(3 * signed_distance_value / epsilon));
for (unsigned int i = 0; i < dim; ++i) {
grad[i] = tau + (1-tau) * factor * grad_signed_distance[i];
}
return grad;
}
private:
double epsilon;
double tau = 1e-8;
SignedDistanceType signed_distance;
};
template <int dim, typename SignedDistanceType>
class NormOfGradient : public Function<dim> {
public:
const PhaseFieldFunction<dim, SignedDistanceType>& pf_function;
NormOfGradient(const PhaseFieldFunction<dim, SignedDistanceType>& pf_function) : Function<dim>(), pf_function(pf_function) {}
virtual double value(const Point<dim> &p, const unsigned int component = 0) const override {
assert(component == 0);
Tensor<1, dim> gradient = pf_function.gradient(p);
return gradient.norm();
}
};
template <int dim, typename SignedDistanceType>
class ExactSolution : public Function<dim> {
public:
ExactSolution(double eps, SignedDistanceType signed_distance)
: Function<dim>(), epsilon(eps), signed_distance(signed_distance), phase_field_function(epsilon, signed_distance) {}
virtual double value(const Point<dim> &p, const unsigned int component = 0) const override {
return phase_field_function.value(p)*(-1/(std::pow(2,dim))*(p.norm_square()-1));
}
Tensor<1, dim> gradient(const Point<dim> &p, const unsigned int component = 0) const override {
Tensor<1, dim> result;
Tensor<1, dim> grad_phase_field = phase_field_function.gradient(p);
double common_term = -1/(std::pow(2,dim))*(p.norm_square()-1);
if (dim == 2) {
result[0] = grad_phase_field[0] * common_term + phase_field_function.value(p) * -1/(std::pow(2,dim)) * 2 * p[0];
result[1] = grad_phase_field[1] * common_term + phase_field_function.value(p) * -1/(std::pow(2,dim)) * 2 * p[1];
} else {
result[0] = grad_phase_field[0] * common_term + phase_field_function.value(p) * -1/(std::pow(2,dim)) * 2 * p[0];
result[1] = grad_phase_field[1] * common_term + phase_field_function.value(p) * -1/(std::pow(2,dim)) * 2 * p[1];
result[2] = grad_phase_field[2] * common_term + phase_field_function.value(p) * -1/(std::pow(2,dim)) * 2 * p[2];
}
return result;
}
private:
SignedDistanceType signed_distance;
PhaseFieldFunction<dim, SignedDistanceType> phase_field_function;
double epsilon;
};
template <int dim, typename SignedDistanceType>
class LaplaceSolver
{
public:
LaplaceSolver(SignedDistanceType signed_distance,
const unsigned int degree,
int h,
const std::vector<double>& eps_list,
double r,
double l,
int lvl_refinement,
const std::string ctrl_mode);
void run();
private:
void make_grid(int cycle);
void refine_mesh();
void setup_discrete_level_set();
void distribute_dofs();
void initialize_matrices();
void assemble_system(int cycle);
double solve();
void output_results(const unsigned int cycle) const;
std::tuple<double, double> process_errors();
bool face_has_ghost_penalty(
const typename Triangulation<dim>::active_cell_iterator &cell,
const unsigned int face_index) const;
std::vector<double> epsilon_list;
int initial_refinement;
unsigned int refinement_levels;
double epsilon;
double radius;
double length;
const std::string weak_formulation;
const unsigned int fe_degree;
const Functions::ConstantFunction<dim> rhs_function;
const Functions::ConstantFunction<dim> boundary_condition;
Triangulation<dim> triangulation;
const FE_Q<dim> fe_level_set;
DoFHandler<dim> level_set_dof_handler;
Vector<double> level_set;
hp::FECollection<dim> fe_collection;
DoFHandler<dim> dof_handler;
Vector<double> solution;
AffineConstraints<double> constraints;
NonMatching::MeshClassifier<dim> mesh_classifier;
SparsityPattern sparsity_pattern;
SparseMatrix<double> stiffness_matrix;
Vector<double> rhs;
SignedDistanceType signed_distance;
PhaseFieldFunction<dim, SignedDistanceType> phase_field;
};
template <int dim, typename SignedDistanceType>
LaplaceSolver<dim, SignedDistanceType>::LaplaceSolver(SignedDistanceType signed_distance,
const unsigned int degree,
int h,
const std::vector<double>& eps_list,
double r,
double l,
int lvl_refinement,
const std::string ctrl_mode)
: fe_degree(degree)
, rhs_function(1.0)
, boundary_condition(0.)
, fe_level_set(fe_degree)
, level_set_dof_handler(triangulation)
, dof_handler(triangulation)
, mesh_classifier(level_set_dof_handler, level_set)
, signed_distance(signed_distance)
, phase_field(eps_list[0], signed_distance)
, initial_refinement(h)
, epsilon_list(eps_list)
, radius(r)
, length(l)
, refinement_levels(lvl_refinement)
, weak_formulation(ctrl_mode)
{}
template <int dim, typename SignedDistanceType>
void LaplaceSolver<dim, SignedDistanceType>::make_grid(int cycle){
std::cout << "Creating background mesh" << std::endl;
triangulation.clear();
GridGenerator::hyper_cube(triangulation, -length/2, length/2);
triangulation.refine_global(initial_refinement+cycle);
}
template <int dim, typename SignedDistanceType>
void LaplaceSolver<dim, SignedDistanceType>::refine_mesh(){
for (unsigned int refinement_cycle = 0; refinement_cycle < refinement_levels; ++refinement_cycle) {
Point<dim> cell_center;
double distance;
double threshold = 3 * epsilon / (refinement_cycle + 1); //3*epsilon , 1.5*epsilon, 1*epsilon
for (const auto &cell : dof_handler.active_cell_iterators()){
cell_center = cell->center();
distance = std::fabs(signed_distance.value(cell_center));
if (distance < threshold) {
cell->set_refine_flag();
}
}
triangulation.execute_coarsening_and_refinement();
}
}
template <int dim, typename SignedDistanceType>
void LaplaceSolver<dim, SignedDistanceType>::setup_discrete_level_set()
{
std::cout << " | Setting up discrete level set function" << std::endl;
level_set_dof_handler.distribute_dofs(fe_level_set);
level_set.reinit(level_set_dof_handler.n_dofs());
VectorTools::interpolate(level_set_dof_handler,
signed_distance,
level_set);
}
enum ActiveFEIndex
{
lagrange = 0,
nothing = 1
};
template <int dim, typename SignedDistanceType>
void LaplaceSolver<dim, SignedDistanceType>::distribute_dofs()
{
fe_collection.push_back(FE_Q<dim>(fe_degree));
fe_collection.push_back(FE_Nothing<dim>());
for (const auto &cell : dof_handler.active_cell_iterators())
{
const NonMatching::LocationToLevelSet cell_location =
mesh_classifier.location_to_level_set(cell);
Point<dim> cell_center = cell->center();
// if (cell_location == NonMatching::LocationToLevelSet::outside)
// if (false)
if (signed_distance.value(cell_center) > 10*epsilon)
cell->set_active_fe_index(ActiveFEIndex::nothing);
else
cell->set_active_fe_index(ActiveFEIndex::lagrange);
}
dof_handler.distribute_dofs(fe_collection);
std::cout << " | Number of degrees of freedom: "
<< dof_handler.n_dofs() << std::endl;
}
template <int dim, typename SignedDistanceType>
void LaplaceSolver<dim, SignedDistanceType>::initialize_matrices()
{
std::cout << " | Initializing matrices" << std::endl;
const auto face_has_flux_coupling = [&](const auto & cell,
const unsigned int face_index) {
return this->face_has_ghost_penalty(cell, face_index);
};
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
const unsigned int n_components = fe_collection.n_components();
Table<2, DoFTools::Coupling> cell_coupling(n_components, n_components);
Table<2, DoFTools::Coupling> face_coupling(n_components, n_components);
cell_coupling[0][0] = DoFTools::always;
face_coupling[0][0] = DoFTools::always;
const bool keep_constrained_dofs = true;
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
constraints.close();
// DoFTools::make_flux_sparsity_pattern(dof_handler,
// dsp,
// constraints,
// keep_constrained_dofs,
// cell_coupling,
// face_coupling,
// numbers::invalid_subdomain_id,
// face_has_flux_coupling);
DoFTools::make_flux_sparsity_pattern(dof_handler,
dsp,
constraints,
keep_constrained_dofs,
numbers::invalid_subdomain_id);
sparsity_pattern.copy_from(dsp);
stiffness_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
rhs.reinit(dof_handler.n_dofs());
}
template <int dim, typename SignedDistanceType>
bool LaplaceSolver<dim, SignedDistanceType>::face_has_ghost_penalty(
const typename Triangulation<dim>::active_cell_iterator &cell,
const unsigned int face_index) const
{
if (cell->at_boundary(face_index))
return false;
Point<dim> cell_center = cell->center();
if (signed_distance.value(cell_center) > 0 && signed_distance.value(cell_center) < 1.5*epsilon){
return true;
}
else return false;
// return false;
}
template <int dim, typename SignedDistanceType>
void LaplaceSolver<dim, SignedDistanceType>::assemble_system(int cycle)
{
std::cout << " | Assembling system" << std::endl;
std::cout << " | Using the " << weak_formulation << " weak formulation" << std::endl;
const unsigned int n_dofs_per_cell = fe_collection[0].dofs_per_cell;
FullMatrix<double> local_stiffness(n_dofs_per_cell, n_dofs_per_cell);
Vector<double> local_rhs(n_dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(n_dofs_per_cell);
const double ghost_parameter = 0.1;
const double nitsche_parameter = 5 * (fe_degree + 1) * fe_degree;
const double h_max = length/std::pow(2,initial_refinement+cycle);
std::cout << " | h_max = " << h_max << std::endl;
const QGauss<dim - 1> face_quadrature(fe_degree + 1);
FEInterfaceValues<dim> fe_interface_values(fe_collection[0],
face_quadrature,
update_gradients |
update_JxW_values |
update_normal_vectors);
const QGauss<dim> quadrature_formula(fe_degree + 1);
FEValues<dim> inside_fe_values(fe_collection[0], quadrature_formula, update_values | update_gradients | update_JxW_values | update_quadrature_points);
double max_gradient_norm = 0;
for (const auto &cell : dof_handler.active_cell_iterators() | IteratorFilters::ActiveFEIndexEqualTo(ActiveFEIndex::lagrange)){
inside_fe_values.reinit(cell);
for (const unsigned int q : inside_fe_values.quadrature_point_indices()){
const Point<dim> &point = inside_fe_values.quadrature_point(q);
if (phase_field.gradient(point).norm() > max_gradient_norm) max_gradient_norm = phase_field.gradient(point).norm();
}
}
for (const auto &cell :
dof_handler.active_cell_iterators() |
IteratorFilters::ActiveFEIndexEqualTo(ActiveFEIndex::lagrange))
{
local_stiffness = 0;
local_rhs = 0;
const double cell_side_length = cell->minimum_vertex_distance();
inside_fe_values.reinit(cell);
if(weak_formulation == "pure_mDDM3"){
for (const unsigned int q : inside_fe_values.quadrature_point_indices()){
const Point<dim> &point = inside_fe_values.quadrature_point(q);
Tensor<1, dim> normal;
if (point.norm() <= radius) normal = - phase_field.gradient(point)/phase_field.gradient(point).norm();
else normal = 0;
for (const unsigned int i : inside_fe_values.dof_indices()){
for (const unsigned int j : inside_fe_values.dof_indices()){
local_stiffness(i, j) +=
//working formulation for mDDM3
(phase_field.value(point) *
inside_fe_values.shape_grad(i, q) *
inside_fe_values.shape_grad(j, q) *
inside_fe_values.JxW(q) +
std::pow(epsilon,-2) *
phase_field.gradient(point).norm() *
inside_fe_values.shape_value(i, q) *
(inside_fe_values.shape_value(j, q) - 0 - signed_distance.value(point) * normal * inside_fe_values.shape_grad(j, q)) *
inside_fe_values.JxW(q));
}
local_rhs(i) +=
// working formulation for mDDM3
(phase_field.value(point) *
1.0 *
inside_fe_values.shape_value(i, q) *
inside_fe_values.JxW(q));
}
}
}
else if (weak_formulation == "Nitsche"){
for (const unsigned int q : inside_fe_values.quadrature_point_indices()){
const Point<dim> &point = inside_fe_values.quadrature_point(q);
Tensor<1, dim> normal;
if (point.norm() <= radius) normal = - phase_field.gradient(point)/phase_field.gradient(point).norm();
else normal = 0;
for (const unsigned int i : inside_fe_values.dof_indices()){
for (const unsigned int j : inside_fe_values.dof_indices()){
local_stiffness(i, j) +=
//Nitsche
(phase_field.value(point) *
inside_fe_values.shape_grad(i, q) *
inside_fe_values.shape_grad(j, q) *
inside_fe_values.JxW(q) -
phase_field.gradient(point).norm()/max_gradient_norm *
inside_fe_values.shape_value(i, q) *
normal * inside_fe_values.shape_grad(j, q) *
inside_fe_values.JxW(q) +
phase_field.gradient(point).norm()/max_gradient_norm *
(std::pow(epsilon,-3) * inside_fe_values.shape_value(i, q) - normal * inside_fe_values.shape_grad(i, q)) *
(inside_fe_values.shape_value(j, q) - 0 - signed_distance.value(point) * normal * inside_fe_values.shape_grad(j, q)) *
inside_fe_values.JxW(q));
}
local_rhs(i) +=
// working formulation for mDDM3
(phase_field.value(point) *
1.0 *
inside_fe_values.shape_value(i, q) *
inside_fe_values.JxW(q));
}
}
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(
local_stiffness, local_rhs, local_dof_indices, stiffness_matrix, rhs);
for (const unsigned int face_index : cell->face_indices())
if (face_has_ghost_penalty(cell, face_index))
{
const unsigned int invalid_subface = numbers::invalid_unsigned_int;
const auto face = cell->face(face_index);
if (face->has_children()){
for (unsigned int subface_no = 0; subface_no < face->n_active_descendants(); ++subface_no){
const auto neighbor_subface_cell = cell->neighbor_child_on_subface(face_index, subface_no);
const auto face_index_neighbor = cell->neighbor_of_neighbor(face_index);
fe_interface_values.reinit(cell,
face_index,
subface_no,
neighbor_subface_cell,
face_index_neighbor,
invalid_subface);
const unsigned int n_interface_dofs =
fe_interface_values.n_current_interface_dofs();
FullMatrix<double> local_stabilization(n_interface_dofs,
n_interface_dofs);
for (unsigned int q = 0; q < fe_interface_values.n_quadrature_points; ++q){
const Tensor<1, dim> normal = fe_interface_values.normal(q);
for (unsigned int i = 0; i < n_interface_dofs; ++i){
for (unsigned int j = 0; j < n_interface_dofs; ++j){
local_stabilization(i, j) +=
.5 * ghost_parameter *
cell_side_length *
normal *
fe_interface_values.jump_in_shape_gradients(i, q) *
normal *
fe_interface_values.jump_in_shape_gradients(j, q) *
fe_interface_values.JxW(q);
}
}
}
const std::vector<types::global_dof_index>
local_interface_dof_indices =
fe_interface_values.get_interface_dof_indices();
constraints.distribute_local_to_global(
local_stabilization, local_interface_dof_indices, stiffness_matrix);
}
}
else{
const auto cell_neighbor = cell->neighbor(face_index);
if (cell->neighbor_is_coarser(face_index)) continue;
if (cell->index() > cell_neighbor->index()){
const auto face_index_neighbor = cell->neighbor_of_neighbor(face_index);
fe_interface_values.reinit(cell,
face_index,
invalid_subface,
cell_neighbor,
face_index_neighbor,
invalid_subface);
const unsigned int n_interface_dofs =
fe_interface_values.n_current_interface_dofs();
FullMatrix<double> local_stabilization(n_interface_dofs,
n_interface_dofs);
for (unsigned int q = 0; q < fe_interface_values.n_quadrature_points; ++q){
const Tensor<1, dim> normal = fe_interface_values.normal(q);
for (unsigned int i = 0; i < n_interface_dofs; ++i){
for (unsigned int j = 0; j < n_interface_dofs; ++j){
local_stabilization(i, j) +=
.5 * ghost_parameter *
cell_side_length *
normal *
fe_interface_values.jump_in_shape_gradients(i, q) *
normal *
fe_interface_values.jump_in_shape_gradients(j, q) *
fe_interface_values.JxW(q);
}
}
}
const std::vector<types::global_dof_index>
local_interface_dof_indices =
fe_interface_values.get_interface_dof_indices();
constraints.distribute_local_to_global(
local_stabilization, local_interface_dof_indices, stiffness_matrix);
}
}
}
}
}
template <int dim, typename SignedDistanceType>
double LaplaceSolver<dim, SignedDistanceType>::solve()
{
std::cout << " | Solving system" << std::endl;
const unsigned int max_iterations = 10000;
SolverControl solver_control(max_iterations, 1e-6);
SolverGMRES<> solver(solver_control);
solver.solve(stiffness_matrix, solution, rhs, PreconditionIdentity());
constraints.distribute(solution);
std::cout << " | | Number of GMRES iterations: " << solver_control.last_step() << std::endl;
std::cout << " | | Final residue of GMRES solver : " << solver_control.last_value() << std::endl;
return solver_control.last_step();
}
template <int dim, typename SignedDistanceType>
void LaplaceSolver<dim, SignedDistanceType>::output_results(const unsigned int cycle) const
{
ExactSolution<dim, SignedDistanceType> exact_solution(epsilon, signed_distance);
NormOfGradient<dim, SignedDistanceType> norm_function(phase_field);
DataOut<dim> data_out;
AffineConstraints<double> projection_constraints;
projection_constraints.clear();
DoFTools::make_hanging_node_constraints(level_set_dof_handler, projection_constraints);
projection_constraints.close();
Vector<double> signed_distance_projection(level_set_dof_handler.n_dofs());
Vector<double> phase_field_projection(level_set_dof_handler.n_dofs());
Vector<double> e_solution(level_set_dof_handler.n_dofs());
Vector<double> has_ghost_penalty(level_set_dof_handler.n_dofs());
Vector<double> gradient_norm(level_set_dof_handler.n_dofs());
VectorTools::project(level_set_dof_handler, projection_constraints, QGauss<dim>(fe_degree+1), signed_distance, signed_distance_projection);
VectorTools::project(level_set_dof_handler, projection_constraints, QGauss<dim>(fe_degree+1), phase_field, phase_field_projection);
VectorTools::project(level_set_dof_handler, projection_constraints, QGauss<dim>(fe_degree+1), exact_solution, e_solution);
VectorTools::project(level_set_dof_handler, projection_constraints, QGauss<dim>(fe_degree+1), norm_function, gradient_norm);
data_out.add_data_vector(dof_handler, solution, "Solution");
data_out.add_data_vector(level_set_dof_handler, signed_distance_projection, "SignedDistance");
data_out.add_data_vector(level_set_dof_handler, phase_field_projection, "PhaseField");
data_out.add_data_vector(level_set_dof_handler, e_solution, "ExactSolution");
data_out.add_data_vector(level_set_dof_handler, gradient_norm, "PhaseFieldGradientNorm");
// data_out.set_cell_selection(
// [this](const typename Triangulation<dim>::cell_iterator &cell) {
// return cell->is_active() && signed_distance.value(cell->center())< epsilon;
// });
data_out.build_patches();
std::string filename = "solution_" + std::to_string(dim) + "D_" + std::to_string(cycle) + ".vtk";
std::ofstream output(filename);
data_out.write_vtk(output);
std::cout << " | | Output written to " << filename << std::endl;
}
template <int dim, typename SignedDistanceType>
std::tuple<double, double> LaplaceSolver<dim, SignedDistanceType>::process_errors(){
ExactSolution<dim, SignedDistanceType> exact_solution(epsilon, signed_distance);
Vector<double> difference_per_cellL2(dof_handler.get_triangulation().n_active_cells());
QGauss<dim> quadrature_formula(fe_level_set.degree + 1); //bilinear quadratic elements
VectorTools::integrate_difference(dof_handler,
solution,
exact_solution,
difference_per_cellL2,
quadrature_formula,
VectorTools::L2_norm); // L2 norm
double L2_error = VectorTools::compute_global_error(dof_handler.get_triangulation(),
difference_per_cellL2,
VectorTools::L2_norm);
Vector<double> difference_per_cellH1(dof_handler.get_triangulation().n_active_cells());
QGauss<dim> quadrature_formulaH1(fe_level_set.degree + 2);
VectorTools::integrate_difference(dof_handler,
solution,
exact_solution,
difference_per_cellH1,
quadrature_formulaH1,
VectorTools::H1_norm); // H1 norm
double H1_error = VectorTools::compute_global_error(dof_handler.get_triangulation(),
difference_per_cellH1,
VectorTools::H1_norm);
return std::make_tuple(L2_error, H1_error);
}
template <int dim, typename SignedDistanceType>
void LaplaceSolver<dim, SignedDistanceType>::run(){
const int num_cycles = epsilon_list.size();
double L2_error, H1_error, convergence_rate_L2, convergence_rate_H1;
std::vector <double> L2_error_list(num_cycles);
std::vector <double> H1_error_list(num_cycles);
std::string output_file;
output_file = "output_" + std::to_string(dim) + "D.csv";
std::ofstream file(output_file);
file << "cycle,epsilon,L2_error,convergence_rate_L2,H1_error,convergence_rate_H1,iterations\n";
double iterations;
for (unsigned int cycle = 0; cycle < num_cycles; ++cycle)
{
std::cout << "\nCycle " << cycle << std::endl;
if (cycle >= num_cycles){
std::cerr << "Error: cycle index exceeds the size of epsilon_list. No output generated for cycle " << std::endl;
return ;
}
else {
epsilon = epsilon_list[cycle];
phase_field.set_epsilon(epsilon);
}
if (cycle == 0)
{
convergence_rate_L2 = 0;
convergence_rate_H1 = 0;
}
make_grid(cycle);
std::cout << " | epsilon = " << epsilon << std::endl;
std::cout << " | Number of active cells: "
<< triangulation.n_active_cells() << std::endl;
refine_mesh();
setup_discrete_level_set();
mesh_classifier.reclassify();
distribute_dofs();
initialize_matrices();
assemble_system(cycle);
iterations = solve();
output_results(cycle);
std::tie(L2_error, H1_error) = process_errors();
L2_error_list[cycle] = L2_error;
H1_error_list[cycle] = H1_error;
if (cycle!=0){
convergence_rate_L2 = log(L2_error_list[cycle]/L2_error_list[cycle-1])/log(epsilon_list[cycle]/epsilon_list[cycle-1]);
convergence_rate_H1 = log(H1_error_list[cycle]/H1_error_list[cycle-1])/log(epsilon_list[cycle]/epsilon_list[cycle-1]);
}
std::cout << " | L2_error = " << L2_error << std::endl;
std::cout << " | convergence_rate_L2 = " << convergence_rate_L2 << std::endl;
std::cout << " | H1_error = " << H1_error << std::endl;
std::cout << " | convergence_rate_H1 = " << convergence_rate_H1 << std::endl;
file << cycle << "," << epsilon << "," << L2_error << "," << convergence_rate_L2 << "," << H1_error << "," << convergence_rate_H1 << "," << iterations << "\n";
}
}
} // namespace Step85
int main()
{
const int dim = 2;
const int degree = 1;
int initial_refinement;
unsigned int refinement_levels = 3;
double radius = 1.;
double length = 2.42;
double h;
const dealii::Functions::SignedDistance::Sphere<dim> signed_distance_sphere;
const std::string weak_formulation = "pure_mDDM3";
// const std::string weak_formulation = "Nitsche";
// const std::vector<double> epsilon_list = {5e-2, 5e-2, 5e-2, 5e-2, 5e-2};
// const std::vector<double> epsilon_list = {2e-1, 1e-1, 5e-2, 2.5e-2};
const std::vector<double> epsilon_list = {2e-1, 1e-1, 5e-2, 2.5e-2, 1.25e-2};
// const std::vector<double> epsilon_list = {2e-1, 1e-1, 5e-2, 2.5e-2, 1.25e-2, 6.25e-3};
// std::vector<double> epsilon_list = {2e-1, 1e-1, 5e-2, 2.5e-2, 1.25e-2, 6.25e-3, 3.13e-3, 1.56e-3};
h = 4.0/16.0;
initial_refinement = static_cast<int>(std::floor(std::log2(length/h)));
try {
std::cout << "\nValues being passed to class instanciation\n -> degree : " << degree << "\n -> initial_refinement : " << initial_refinement << "\n -> radius : " << radius << std::endl;
std::cout << " -> length : " << length << "\n -> refinement_levels : " << refinement_levels << std::endl;
Step85::LaplaceSolver<dim, dealii::Functions::SignedDistance::Sphere<dim>> laplace_solver(signed_distance_sphere,
degree,
initial_refinement,
epsilon_list, radius,
length,
refinement_levels,
weak_formulation);
laplace_solver.run();
std::cout << "\n";
}
catch (std::exception &exc){
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...){
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
}
