Dear Jean-Paul, Thank you very much for your reply. Please find the sample code that reproduce the problem. The discontinuity and overlap can be seen by zooming into the mid-bottom and mid-upper part of the deformed configuration.
Thanks, Anup. On Saturday, January 21, 2017 at 4:53:13 AM UTC-6, Jean-Paul Pelteret wrote: > > Dear Anup, > > Its possible that this is related to this issue > <https://github.com/dealii/dealii/issues/2206>. Do you perhaps have a > minimal test case that reproduces the problem? > > Best, > Jean-Paul > > On Saturday, January 21, 2017 at 4:07:00 AM UTC+1, Anup Basak wrote: >> >> Hello all, >> >> I have solved a 2D mechanics (thermoelasticity) problem based on step 44 >> using quadratic elements (FE_Q). When I plot the solution in the deformed >> configuration (attached herewith) I see that at one place discontinuity >> gets developed between two elements (Fig. in page 2) and at some other >> place two elements overlap (Fig. in page 3). I have shown the locations in >> the figure of page 1. However, when linear elements are used, this problem >> does not appear. In both the cases I get sufficient convergence. I have >> confirmed this issue with both Visit and Paraview. I never encountered such >> strange problem while solving 3D problems. Hence it might be an issue >> related to post-processing. >> >> Has anyone encountered such problem in the past? I shall be thankful if >> someone could tell me what might go wrong. >> >> Thanks and regards, >> Anup. >> > -- 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 dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
#include <deal.II/base/function.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/base/std_cxx11/shared_ptr.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/fe/fe_dgp_monomial.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q_eulerian.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition_selector.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <iostream>
#include <fstream>
//////////////////////////////////////////////////////////////////
namespace PhaseField
{
using namespace dealii;
// INPUT OF PARAMETERS
namespace Parameters
{
struct FESystem
{
unsigned int poly_degree;
unsigned int quad_order;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void FESystem::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Finite element system");
{
prm.declare_entry("Polynomial degree", "2",
Patterns::Integer(0),
"Displacement system polynomial order");
prm.declare_entry("Quadrature order", "3",
Patterns::Integer(0),
"Gauss quadrature order");
}
prm.leave_subsection();
}
void FESystem::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Finite element system");
{
poly_degree = prm.get_integer("Polynomial degree");
quad_order = prm.get_integer("Quadrature order");
}
prm.leave_subsection();
}
////////////////////////////////////////////////////
struct Geometry
{
unsigned int elements_per_width;
double width;
double scale;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void Geometry::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Geometry");
{
prm.declare_entry("Elements per width", "40",
Patterns::Integer(0),
"Number of elements per width");
prm.declare_entry("width", "50.0",
Patterns::Double(0),
"width");
prm.declare_entry("Grid scale", "1e-9",
Patterns::Double(0.0),
"Global grid scaling factor");
}
prm.leave_subsection();
}
void Geometry::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Geometry");
{
elements_per_width = prm.get_integer("Elements per width");
width = prm.get_double("width");
scale = prm.get_double("Grid scale");
}
prm.leave_subsection();
}
/////////////////////////////////////////////////
struct Materials
{
double lambda0;
double mu0;
double lambdai;
double mui;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void Materials::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Material properties");
{
prm.declare_entry("lambda austenite", "0.0", /* */
Patterns::Double(),
"lambda austenite");
prm.declare_entry("mu austenite", "0.0",
Patterns::Double(0.0),
"mu austenite");
prm.declare_entry("lambda martensite", "0.0",
Patterns::Double(),
"lambda martensite");
prm.declare_entry("mu martensite", "0.0",
Patterns::Double(0.0),
"mu martensite");
}
prm.leave_subsection();
}
void Materials::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Material properties");
{
lambda0 = prm.get_double("lambda austenite");
mu0 = prm.get_double("mu austenite");
lambdai = prm.get_double("lambda martensite");
mui = prm.get_double("mu martensite");
}
prm.leave_subsection();
}
////////////////////////////////////////////////////////
struct LinearSolver
{
std::string type_lin;
double tol_lin;
double max_iterations_lin;
std::string preconditioner_type;
double preconditioner_relaxation;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void LinearSolver::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Linear solver");
{
prm.declare_entry("Solver type", "CG",
Patterns::Selection("CG|Direct"),
"Type of solver used to solve the linear system");
prm.declare_entry("Residual", "1e-12",
Patterns::Double(0.0),
"Linear solver residual (scaled by residual norm)");
prm.declare_entry("Max iteration multiplier", "1",
Patterns::Double(0.0),
"Linear solver iterations (multiples of the system matrix size)");
prm.declare_entry("Preconditioner type", "ssor",
Patterns::Selection("jacobi|ssor"),
"Type of preconditioner");
prm.declare_entry("Preconditioner relaxation", "0.65",
Patterns::Double(0.0),
"Preconditioner relaxation value");
}
prm.leave_subsection();
}
void LinearSolver::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Linear solver");
{
type_lin = prm.get("Solver type");
tol_lin = prm.get_double("Residual");
max_iterations_lin = prm.get_double("Max iteration multiplier");
preconditioner_type = prm.get("Preconditioner type");
preconditioner_relaxation = prm.get_double("Preconditioner relaxation");
}
prm.leave_subsection();
}
///////////////////////////////////////////////////
struct NonlinearSolver
{
unsigned int max_iterations_NR;
double tol_f;
double tol_u;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void NonlinearSolver::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Nonlinear solver");
{
prm.declare_entry("Max iterations Newton-Raphson", "10",
Patterns::Integer(0),
"Number of Newton-Raphson iterations allowed");
prm.declare_entry("Tolerance force", "1.0e-12",
Patterns::Double(0.0),
"Force residual tolerance");
prm.declare_entry("Tolerance displacement", "1.0e-12",
Patterns::Double(0.0),
"Displacement error tolerance");
}
prm.leave_subsection();
}
void NonlinearSolver::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Nonlinear solver");
{
max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
tol_f = prm.get_double("Tolerance force");
tol_u = prm.get_double("Tolerance displacement");
}
prm.leave_subsection();
}
/////////////////////////////////////////////////
struct Time
{
double delta_t;
double end_time;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void Time::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Time");
{
prm.declare_entry("End time", "0.0",
Patterns::Double(),
"End time");
prm.declare_entry("Time step size", "0.0",
Patterns::Double(),
"Time step size");
}
prm.leave_subsection();
}
void Time::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Time");
{
end_time = prm.get_double("End time");
delta_t = prm.get_double("Time step size");
}
prm.leave_subsection();
}
///////////////////////////////////////////////////////
struct AllParameters : public FESystem,
public Geometry,
public Materials,
public LinearSolver,
public NonlinearSolver,
public Time
{
AllParameters(const std::string &input_file);
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
AllParameters::AllParameters(const std::string &input_file)
{
ParameterHandler prm;
declare_parameters(prm);
prm.read_input(input_file);
parse_parameters(prm);
}
void AllParameters::declare_parameters(ParameterHandler &prm)
{
FESystem::declare_parameters(prm);
Geometry::declare_parameters(prm);
Materials::declare_parameters(prm);
LinearSolver::declare_parameters(prm);
NonlinearSolver::declare_parameters(prm);
Time::declare_parameters(prm);
}
void AllParameters::parse_parameters(ParameterHandler &prm)
{
FESystem::parse_parameters(prm);
Geometry::parse_parameters(prm);
Materials::parse_parameters(prm);
LinearSolver::parse_parameters(prm);
NonlinearSolver::parse_parameters(prm);
Time::parse_parameters(prm);
}
}
// DEFINE SECOND ORDER IDENTITY, AND TWO FOURTH ORDER IDENTITY TENSORS
template <int dim>
class StandardTensors
{
public:
static const SymmetricTensor<2, dim> I;
static const SymmetricTensor<4, dim> IxI;
static const SymmetricTensor<4, dim> II;
};
template <int dim>
const SymmetricTensor<2, dim>
StandardTensors<dim>::I = unit_symmetric_tensor<dim>();
template <int dim>
const SymmetricTensor<4, dim>
StandardTensors<dim>::IxI = outer_product(I, I);
template <int dim>
const SymmetricTensor<4, dim>
StandardTensors<dim>::II = identity_tensor<dim>();
// DEFINE TIME STEP, CURRENT TIME ETC.
class Time
{
public:
Time (const double time_end,
const double delta_t)
:
timestep(0),
time_current(0.0),
time_end(time_end),
delta_t(delta_t),
delta_t_adapt(delta_t)
{}
virtual ~Time()
{}
double current() const
{
return time_current;
}
double end() const
{
return time_end;
}
double get_delta_t() const
{
return delta_t_adapt;
}
unsigned int get_timestep() const
{
return timestep;
}
void increment()
{
time_current += delta_t_adapt;
++timestep;
}
private:
unsigned int timestep;
double time_current;
const double time_end;
const double delta_t;
double delta_t_adapt;
};
//////////////////////////////////////////////////////////
template <int dim>
class Material_Compressible_Neo_Hook_Three_Field
{
public:
Material_Compressible_Neo_Hook_Three_Field(const double lambda0,
const double mu0,
const double lambdai,
const double mui)
:
det_F(1.0),
I1(0.0),
ge(StandardTensors<dim>::I),
Ee(SymmetricTensor<2, dim>()),
lambda0(lambda0),
mu0(mu0),
lambdai(lambdai),
mui(mui)
{}
~Material_Compressible_Neo_Hook_Three_Field()
{}
void update_material_data (const Tensor<2, dim> &F, const Tensor<2, dim> &Fe,
const double det_Ft, const double Ut_22, const double etai, const double phii)
{
det_F = determinant(F);
detFt = det_Ft;
Ce = symmetrize(transpose(Fe)*Fe);
ge = symmetrize(Fe*transpose(Fe));
Ee = 0.5*(Ce-StandardTensors<dim>::I);
I1 = trace(Ee);
lambda00 = 2.0*mu0*lambda0/(lambda0+2.0*mu0);
lambdaii = 2.0*mui*lambdai/(lambdai+2.0*mui);
lambda = lambda00;
mu = mu0;
double Fe_22;
double Fe_22_p = sqrt(1.0-lambda/mu*I1);
if (det_F*Fe_22_p*Ut_22 > 0.)
Fe_22 = Fe_22_p;
else
Fe_22 = -Fe_22_p;
F_22 = Fe_22*Ut_22;
Assert(det_F > 0, ExcInternalError());
}
SymmetricTensor<2, dim> get_tau()
{
return get_tau_kirchhoff();
}
SymmetricTensor<4, dim> get_Jc() const
{
return get_Jc_modulus();
}
double get_det_F() const
{
return det_F;
}
protected:
double F_22;
double det_F;
double detFt;
double I1;
SymmetricTensor<2, dim> ge;
SymmetricTensor<2, dim> Ce;
SymmetricTensor<2, dim> Ee;
double lambda0;
double mu0;
double lambdai;
double mui;
double lambda00;
double mu00;
double lambdaii;
double muii;
double lambda;
double mu;
// compute the Kirchhoff stress
SymmetricTensor<2, dim> get_tau_kirchhoff() const
{
SymmetricTensor<2, dim> kirchhoff_stress =
(lambda*I1*ge + mu*symmetrize(Tensor<2, dim>(ge) * Tensor<2, dim>(ge))-mu*ge)*detFt;
return kirchhoff_stress;
}
// compute the modulus J*d_ijkl
SymmetricTensor<4, dim> get_Jc_modulus() const
{
SymmetricTensor<4,dim> elasticityTensor;
for (unsigned int i=0; i<dim; ++i)
for (unsigned int j=0; j<dim; ++j)
for (unsigned int k=0; k<dim; ++k)
for (unsigned int l=0; l<dim; ++l)
elasticityTensor[i][j][k][l]
= (lambda*ge[i][j]*ge[k][l]+mu*(ge[i][l]*ge[j][k]+ge[i][k]*ge[j][l]));
return elasticityTensor;
}
};
///////////////////////// updates the quadrature point history
template <int dim>
class PointHistory
{
public:
PointHistory()
:
material(NULL),
F_inv(StandardTensors<dim>::I),
tau(SymmetricTensor<2, dim>()),
Jc(SymmetricTensor<4, dim>())
{}
virtual ~PointHistory()
{
delete material;
material = NULL;
}
void setup_lqp (const Parameters::AllParameters ¶meters)
{
material = new Material_Compressible_Neo_Hook_Three_Field<dim>(parameters.lambda0,
parameters.mu0, parameters.lambdai, parameters.mui);
update_values(Tensor<2, dim>(), double ());
}
void update_values (const Tensor<2, dim> &Grad_u_n, const double etai)
{
const Tensor<2, dim> F = (Tensor<2, dim>(StandardTensors<dim>::I) + Grad_u_n);
const double phii = 3.0*etai*etai-2.0*etai*etai*etai;
Tensor<2, dim> et1;
Tensor<2, dim> et2;
et1[0][0] = 0.0285;
et1[0][1] = -0.0665;
et1[1][0] = et1[0][1];
et1[1][1] = 0.0285;
et2[0][0] = 0.0285;
et2[0][1] = 0.0665;
et2[1][0] = et2[0][1];
et2[1][1] = 0.0285;
Ft = Tensor<2, dim>(StandardTensors<dim>::I) + et2 + (et1-et2)*phii;
double Ut_22 = 0.922;
det_Ft = determinant(Ft)*Ut_22;
Fe = F * invert(Ft);
material->update_material_data(F, Fe, det_Ft, Ut_22, etai, phii);
F_inv = invert(F);
tau = material->get_tau();
Jc = material->get_Jc();
Assert(determinant(F_inv) > 0, ExcInternalError());
}
double get_det_F() const
{
return material->get_det_F();
}
double get_det_Ft() const
{
return det_Ft;
}
const Tensor<2, dim> &get_F_inv() const
{
return F_inv;
}
const Tensor<2, dim> &get_Ft() const
{
return Ft;
}
const SymmetricTensor<2, dim> &get_tau() const
{
return tau;
}
const SymmetricTensor<4, dim> &get_Jc() const
{
return Jc;
}
double get_phi_tilde() const
{
return phi_tilde;
}
private:
Material_Compressible_Neo_Hook_Three_Field<dim> *material;
Tensor<2, dim> F_inv;
Tensor<2, dim> Ft;
SymmetricTensor<2, dim> tau;
SymmetricTensor<4, dim> Jc;
Tensor<2, dim> Fe;
double det_Ft;
double phi_tilde;
};
///////////////////////////////////////////////////////////////
template <int dim>
class Solid
{
public:
Solid(const std::string &input_file);
virtual
~Solid();
void
run();
private:
struct PerTaskData_K;
struct ScratchData_K;
struct PerTaskData_RHS;
struct ScratchData_RHS;
struct PerTaskData_UQPH;
struct ScratchData_UQPH;
void
make_grid();
void
system_setup();
void
assemble_system_tangent();
void
assemble_system_tangent_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_K &scratch,
PerTaskData_K &data);
void
copy_local_to_global_K(const PerTaskData_K &data);
void
assemble_system_rhs();
void
assemble_system_rhs_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_RHS &scratch,
PerTaskData_RHS &data);
void
copy_local_to_global_rhs(const PerTaskData_RHS &data);
void
make_constraints(const int &it_nr);
void
setup_qph();
void update_qph_incremental(const BlockVector<double> &solution_delta);
void
update_qph_incremental_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_UQPH &scratch,
PerTaskData_UQPH &data);
void
copy_local_to_global_UQPH(const PerTaskData_UQPH &/*data*/)
{}
void
solve_nonlinear_timestep(BlockVector<double> &solution_delta);
// performs the Newton-Raphson iteration to compute the incremental displacements
std::pair<unsigned int, double>
solve_linear_system(BlockVector<double> &newton_update);
// computes total displacements
BlockVector<double>
get_total_solution(const BlockVector<double> &solution_delta) const;
// assemble matrices related to eta
void
output_results() const;
Parameters::AllParameters parameters;
double vol_reference; // volume of the reference config
double vol_current; // volume of the current config
Triangulation<dim> triangulation;
Time time; // variable of type class 'Time'
TimerOutput timer;
double delta_t_adapt;
std::vector<PointHistory<dim> > quadrature_point_history;
const unsigned int degree; // degree of polynomial of shape functions
const FESystem<dim> fe; // fe object
DoFHandler<dim> dof_handler_ref;
const unsigned int dofs_per_cell; // no of dofs per cell for the mechanics problem
const FEValuesExtractors::Vector u_fe; // this has been used from step-44, this can be avoided
// if we eliminate the block structure
static const unsigned int n_blocks = 1; // no of block structure
static const unsigned int n_components = dim; // no of displacement components
static const unsigned int first_u_component = 0;
enum
{
u_dof = 0
};
std::vector<types::global_dof_index> dofs_per_block; // total no of dofs per block
const QGauss<dim> qf_cell; // quadrature points in the cell
const QGauss<dim - 1> qf_face; // quadrature points at the face
const unsigned int n_q_points; // no of quadrature points in the cell
const unsigned int n_q_points_f; // no of quadrature points at the face
ConstraintMatrix constraints; // constraint object
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> tangent_matrix; // tangent stiffenss matrix
BlockVector<double> system_rhs; // system right hand side or residual of mechanics problem
BlockVector<double> solution_n; // solution vector for displacement
BlockVector<double> solution_n_iter; // another vector containing the displacement soln
const unsigned int degree_eta; // degree of polynomial for eta
FE_Q<dim> fe_eta; // fe object for eta
DoFHandler<dim> dof_handler_eta; //another dof_handler for eta
const unsigned int dofs_per_cell_eta; // dof per eta cell
SparsityPattern sparsity_pattern_eta;
Vector<double> solution_etai; // solution vector for eta
// structure defining various errors
struct Errors
{
Errors()
:
norm(1.0), u(1.0)
{}
void reset()
{
norm = 1.0;
u = 1.0;
}
void normalise(const Errors &rhs)
{
if (rhs.norm != 0.0)
norm /= rhs.norm;
if (rhs.u != 0.0)
u /= rhs.u;
}
double norm, u;
};
Errors error_residual, error_residual_0, error_residual_norm, error_update,
error_update_0, error_update_norm;
void
get_error_residual(Errors &error_residual);
void
get_error_update(const BlockVector<double> &newton_update,
Errors &error_update);
static
void
print_conv_header();
void
print_conv_footer();
};
/////////////////////////////////////////////////////////
// defines the initial condition for the order parameter
template <int dim>
class Initial_etai : public Function<dim>
{
public:
Initial_etai () : Function<dim>() {}
virtual double value(const Point<dim> &p,
const unsigned int /*component = 0*/) const;
};
template <int dim>
double Initial_etai<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
if (p[0] < 10.0e-9)
return 1.0;
else if (p[0] > 10.0e-9)
return 0.0;
}
///////////////////////////////////////////////////
// constructor
template <int dim>
Solid<dim>::Solid(const std::string &input_file)
:
parameters(input_file),
triangulation(Triangulation<dim>::maximum_smoothing),
time(parameters.end_time, parameters.delta_t),
timer(std::cout,
TimerOutput::summary,
TimerOutput::wall_times),
degree(parameters.poly_degree),
fe(FE_Q<dim>(parameters.poly_degree), dim), // displacement
dof_handler_ref(triangulation),
dofs_per_cell (fe.dofs_per_cell),
u_fe(first_u_component),
dofs_per_block(n_blocks),
qf_cell(parameters.quad_order),
qf_face(parameters.quad_order),
n_q_points (qf_cell.size()),
n_q_points_f (qf_face.size()),
degree_eta(parameters.poly_degree),
fe_eta (parameters.poly_degree),
dof_handler_eta (triangulation),
dofs_per_cell_eta(fe_eta.dofs_per_cell)
{ }
template <int dim>
Solid<dim>::~Solid()
{
dof_handler_ref.clear();
dof_handler_eta.clear();
}
//////////////////////////////////////////////////
template <int dim>
void Solid<dim>::run()
{
make_grid(); // generates the geometry and mesh
system_setup(); // sets up the system matrices
VectorTools::interpolate(dof_handler_eta, Initial_etai<dim>(), solution_etai);//initial eta
output_results(); // prints output
BlockVector<double> solution_delta(dofs_per_block); // du (displacement vector)
solution_delta = 0.0;
solution_n_iter = solution_n;
solve_nonlinear_timestep(solution_delta);
solution_n += solution_delta;
output_results();
}
// Next three functions are used in computing and assembling the system tangent stiffness matrix
/////////////////////////////////////////////////////////////
template <int dim>
struct Solid<dim>::PerTaskData_K
{
FullMatrix<double> cell_matrix;
std::vector<types::global_dof_index> local_dof_indices;
PerTaskData_K(const unsigned int dofs_per_cell)
:
cell_matrix(dofs_per_cell, dofs_per_cell),
local_dof_indices(dofs_per_cell)
{}
void reset()
{
cell_matrix = 0.0;
}
};
/////////////////////////////////////////////////////////
template <int dim>
struct Solid<dim>::ScratchData_K
{
FEValues<dim> fe_values_ref;
std::vector<std::vector<double> > Nx;
std::vector<std::vector<Tensor<2, dim> > > grad_Nx;
std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
ScratchData_K(const FiniteElement<dim> &fe_cell,
const QGauss<dim> &qf_cell,
const UpdateFlags uf_cell)
:
fe_values_ref(fe_cell, qf_cell, uf_cell),
Nx(qf_cell.size(),
std::vector<double>(fe_cell.dofs_per_cell)),
grad_Nx(qf_cell.size(),
std::vector<Tensor<2, dim> >(fe_cell.dofs_per_cell)),
symm_grad_Nx(qf_cell.size(),
std::vector<SymmetricTensor<2, dim> >
(fe_cell.dofs_per_cell))
{}
ScratchData_K(const ScratchData_K &rhs)
:
fe_values_ref(rhs.fe_values_ref.get_fe(),
rhs.fe_values_ref.get_quadrature(),
rhs.fe_values_ref.get_update_flags()),
Nx(rhs.Nx),
grad_Nx(rhs.grad_Nx),
symm_grad_Nx(rhs.symm_grad_Nx)
{}
void reset()
{
const unsigned int n_q_points = Nx.size();
const unsigned int n_dofs_per_cell = Nx[0].size();
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
Assert( grad_Nx[q_point].size() == n_dofs_per_cell,
ExcInternalError());
Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell,
ExcInternalError());
for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
{
Nx[q_point][k] = 0.0;
grad_Nx[q_point][k] = 0.0;
symm_grad_Nx[q_point][k] = 0.0;
}
}
}
};
//////////////////////////////////////////////////////////////
template <int dim>
struct Solid<dim>::PerTaskData_RHS
{
Vector<double> cell_rhs;
std::vector<types::global_dof_index> local_dof_indices;
PerTaskData_RHS(const unsigned int dofs_per_cell)
:
cell_rhs(dofs_per_cell),
local_dof_indices(dofs_per_cell)
{}
void reset()
{
cell_rhs = 0.0;
}
};
template <int dim>
struct Solid<dim>::ScratchData_RHS
{
FEValues<dim> fe_values_ref;
FEFaceValues<dim> fe_face_values_ref;
std::vector<std::vector<double> > Nx;
std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
ScratchData_RHS(const FiniteElement<dim> &fe_cell,
const QGauss<dim> &qf_cell, const UpdateFlags uf_cell,
const QGauss<dim - 1> &qf_face, const UpdateFlags uf_face)
:
fe_values_ref(fe_cell, qf_cell, uf_cell),
fe_face_values_ref(fe_cell, qf_face, uf_face),
Nx(qf_cell.size(),
std::vector<double>(fe_cell.dofs_per_cell)),
symm_grad_Nx(qf_cell.size(),
std::vector<SymmetricTensor<2, dim> >
(fe_cell.dofs_per_cell))
{}
ScratchData_RHS(const ScratchData_RHS &rhs)
:
fe_values_ref(rhs.fe_values_ref.get_fe(),
rhs.fe_values_ref.get_quadrature(),
rhs.fe_values_ref.get_update_flags()),
fe_face_values_ref(rhs.fe_face_values_ref.get_fe(),
rhs.fe_face_values_ref.get_quadrature(),
rhs.fe_face_values_ref.get_update_flags()),
Nx(rhs.Nx),
symm_grad_Nx(rhs.symm_grad_Nx)
{}
void reset()
{
const unsigned int n_q_points = Nx.size();
const unsigned int n_dofs_per_cell = Nx[0].size();
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell,
ExcInternalError());
for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
{
Nx[q_point][k] = 0.0;
symm_grad_Nx[q_point][k] = 0.0;
}
}
}
};
/////////////////////////////////////////////////////
template <int dim>
struct Solid<dim>::PerTaskData_UQPH
{
void reset()
{}
};
////////////////////////////////////////////////////
template <int dim>
struct Solid<dim>::ScratchData_UQPH
{
const BlockVector<double> &solution_total;
std::vector<Tensor<2, dim> > solution_grads_u_total;
std::vector<double> solution_etai_total;
FEValues<dim> fe_values_ref;
FEValues<dim> fe_values_eta;
ScratchData_UQPH(const FiniteElement<dim> &fe_cell,
const UpdateFlags uf_cell,
const FiniteElement<dim> &fe_cell_eta,
const UpdateFlags uf_cell_eta,
const QGauss<dim> &qf_cell,
const BlockVector<double> &solution_total)
:
solution_total(solution_total),
solution_grads_u_total(qf_cell.size()),
solution_etai_total (qf_cell.size()),
fe_values_ref(fe_cell, qf_cell, uf_cell),
fe_values_eta(fe_cell_eta, qf_cell, uf_cell_eta)
{}
ScratchData_UQPH(const ScratchData_UQPH &rhs)
:
solution_total(rhs.solution_total),
solution_grads_u_total(rhs.solution_grads_u_total),
solution_etai_total(rhs.solution_etai_total),
fe_values_ref(rhs.fe_values_ref.get_fe(),
rhs.fe_values_ref.get_quadrature(),
rhs.fe_values_ref.get_update_flags()),
fe_values_eta(rhs.fe_values_eta.get_fe(),
rhs.fe_values_eta.get_quadrature(),
rhs.fe_values_eta.get_update_flags())
{}
void reset()
{
const unsigned int n_q_points = solution_grads_u_total.size();
for (unsigned int q = 0; q < n_q_points; ++q)
{
solution_grads_u_total[q] = 0.0;
solution_etai_total[q] = 0.0;
}
}
};
///////////////////////////////////////////////////////
template <int dim>
Point<dim> grid_y_transform (const Point<dim> &pt_in)
{
const double &x = pt_in[0];
Point<dim> pt_out = pt_in;
return pt_out;
}
////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::make_grid()
{
std::vector< unsigned int > repetitions(dim, parameters.elements_per_width);
repetitions[0] = parameters.elements_per_width;
const Point<dim> bottom_left = (dim == 3 ? Point<dim>(0.0, 0.0, -0.5) : Point<dim>(0.0, 0.0));
const Point<dim> top_right = (dim == 3 ? Point<dim>(parameters.width, parameters.width,
0.5) :
Point<dim>(parameters.width, parameters.width));
GridGenerator::subdivided_hyper_rectangle(triangulation,
repetitions,
bottom_left,
top_right);
GridTools::transform(&grid_y_transform<dim>, triangulation);
const double tol_boundary = 1e-6;
typename Triangulation<dim>::active_cell_iterator cell =
triangulation.begin_active(), endc = triangulation.end();
for (; cell != endc; ++cell)
for (unsigned int face = 0;
face < GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary() == true)
{
if (std::abs(cell->face(face)->center()[0] - 0.0) < tol_boundary)
cell->face(face)->set_boundary_id(1); // -X faces
}
GridTools::scale(parameters.scale, triangulation);
vol_reference = GridTools::volume(triangulation);
vol_current = vol_reference;
std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
}
///////////////////////////////////////////////////
template <int dim>
void Solid<dim>::system_setup()
{ //std::cout << " setup " << std::flush;
timer.enter_subsection("Setup system");
std::vector<unsigned int> block_component(n_components, u_dof);
dof_handler_ref.distribute_dofs(fe);
DoFRenumbering::Cuthill_McKee(dof_handler_ref);
DoFRenumbering::component_wise(dof_handler_ref, block_component);
DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block,
block_component);
std::cout << "Triangulation:"
<< "\n\t Number of active cells: " << triangulation.n_active_cells()
<< "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs()
<< std::endl;
// block system for the mechanics problem, but block system can be avoided
tangent_matrix.clear();
{
const types::global_dof_index n_dofs_u = dofs_per_block[u_dof];
BlockDynamicSparsityPattern csp(n_blocks, n_blocks);
csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
csp.collect_sizes();
Table<2, DoFTools::Coupling> coupling(n_components, n_components);
for (unsigned int ii = 0; ii < n_components; ++ii)
for (unsigned int jj = 0; jj < n_components; ++jj)
coupling[ii][jj] = DoFTools::always;
DoFTools::make_sparsity_pattern(dof_handler_ref,
coupling,
csp,
constraints,
false);
sparsity_pattern.copy_from(csp);
}
// initialization of matrices and vectors for the mechanics problem
tangent_matrix.reinit(sparsity_pattern);
system_rhs.reinit(dofs_per_block);
system_rhs.collect_sizes();
solution_n.reinit(dofs_per_block);
solution_n.collect_sizes();
solution_n_iter.reinit(dofs_per_block);
solution_n_iter.collect_sizes();
// initialization of matrices and vectors for the G-L equation
dof_handler_eta.distribute_dofs (fe_eta);
const unsigned int n_dofs_eta = dof_handler_eta.n_dofs();
{
DynamicSparsityPattern dsp (n_dofs_eta, n_dofs_eta);
DoFTools::make_sparsity_pattern (dof_handler_eta, dsp);
sparsity_pattern_eta.copy_from(dsp);
solution_etai.reinit(dof_handler_eta.n_dofs());
}
// sets up the quadrature point history
setup_qph();
timer.leave_subsection();
}
//////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::setup_qph()
{
std::cout << " Setting up quadrature point data..." << std::endl;
{
triangulation.clear_user_data();
{
std::vector<PointHistory<dim> > tmp;
tmp.swap(quadrature_point_history);
}
quadrature_point_history
.resize(triangulation.n_active_cells() * n_q_points);
unsigned int history_index = 0;
for (typename Triangulation<dim>::active_cell_iterator cell =
triangulation.begin_active(); cell != triangulation.end();
++cell)
{
cell->set_user_pointer(&quadrature_point_history[history_index]);
history_index += n_q_points;
}
Assert(history_index == quadrature_point_history.size(),
ExcInternalError());
}
for (typename Triangulation<dim>::active_cell_iterator cell =
triangulation.begin_active(); cell != triangulation.end(); ++cell)
{
PointHistory<dim> *lqph =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
Assert(lqph <= &quadrature_point_history.back(), ExcInternalError());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
lqph[q_point].setup_lqp(parameters);
}
}
// Updates the data at all quadrature points over a loop run by WorkStream::run
template <int dim>
void Solid<dim>::update_qph_incremental(const BlockVector<double> &solution_delta)
{
timer.enter_subsection("Update QPH data");
std::cout << " UQPH " << std::flush;
const BlockVector<double> solution_total(get_total_solution(solution_delta));
const UpdateFlags uf_UQPH(update_gradients);
const UpdateFlags uf_UQPH_eta(update_values | update_gradients);
PerTaskData_UQPH per_task_data_UQPH;
ScratchData_UQPH scratch_data_UQPH(fe, uf_UQPH, fe_eta, uf_UQPH_eta, qf_cell, solution_total);
WorkStream::run(dof_handler_ref.begin_active(),
dof_handler_ref.end(),
*this,
&Solid::update_qph_incremental_one_cell,
&Solid::copy_local_to_global_UQPH,
scratch_data_UQPH,
per_task_data_UQPH);
timer.leave_subsection();
}
// updates the data at the quadrature points in a single cell
template <int dim>
void
Solid<dim>::update_qph_incremental_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_UQPH &scratch,
PerTaskData_UQPH &/*data*/)
{
PointHistory<dim> *lqph =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
Assert(lqph <= &quadrature_point_history.back(), ExcInternalError());
Assert(scratch.solution_grads_u_total.size() == n_q_points,
ExcInternalError());
Assert(scratch.solution_etai_total.size() == n_q_points,
ExcInternalError());
scratch.reset();
scratch.fe_values_ref.reinit(cell);
typename DoFHandler<dim>::active_cell_iterator
eta_cell (&triangulation,
cell->level(),
cell->index(),
&dof_handler_eta);
scratch.fe_values_eta.reinit (eta_cell);
scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total,
scratch.solution_grads_u_total);
scratch.fe_values_eta.get_function_values(solution_etai,
scratch.solution_etai_total);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
lqph[q_point].update_values(scratch.solution_grads_u_total[q_point],
scratch.solution_etai_total[q_point]);
}
// solves the N-R problem for displacement
template <int dim>
void
Solid<dim>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
{
std::cout << std::endl << "Timestep " << time.get_timestep() << " @ "
<< time.current() << "s" << std::endl;
BlockVector<double> newton_update(dofs_per_block);
error_residual.reset();
error_residual_0.reset();
error_residual_norm.reset();
error_update.reset();
error_update_0.reset();
error_update_norm.reset();
print_conv_header();
unsigned int newton_iteration = 0;
for (; newton_iteration < parameters.max_iterations_NR;
++newton_iteration)
{
std::cout << " " << std::setw(2) << newton_iteration << " " << std::flush;
tangent_matrix = 0.0;
system_rhs = 0.0;
assemble_system_rhs(); // assembling the residual
get_error_residual(error_residual);
if (newton_iteration == 0)
error_residual_0 = error_residual;
error_residual_norm = error_residual;
error_residual_norm.normalise(error_residual_0);
// if the norm of residual less than the tolerance, break
if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u
&& error_residual_norm.u <= parameters.tol_f)
{
std::cout << " CONVERGED! " << std::endl;
print_conv_footer();
break;
}
assemble_system_tangent(); // assembles the tangent stiffness
make_constraints(newton_iteration); // applying the constraints
constraints.condense(tangent_matrix, system_rhs);
const std::pair<unsigned int, double>
lin_solver_output = solve_linear_system(newton_update); // newton update is du
get_error_update(newton_update, error_update);
if (newton_iteration == 0)
error_update_0 = error_update;
error_update_norm = error_update;
error_update_norm.normalise(error_update_0);
solution_delta += newton_update;
update_qph_incremental(solution_delta);
std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
<< std::scientific << lin_solver_output.first << " "
<< lin_solver_output.second << " " << error_residual_norm.norm
<< " " << error_residual_norm.u << " "
<< " " << error_update_norm.norm << " " << error_update_norm.u
<< " " << std::endl;
}
AssertThrow (newton_iteration <= parameters.max_iterations_NR,
ExcMessage("No convergence in nonlinear solver!"));
}
///////////////////////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::print_conv_header()
{
static const unsigned int l_width = 102;
for (unsigned int i = 0; i < l_width; ++i)
std::cout << "_";
std::cout << std::endl;
std::cout << " SOLVER STEP "
<< " | LIN_IT LIN_RES RES_NORM "
<< " RES_U NU_NORM "
<< " NU_U " << std::endl;
for (unsigned int i = 0; i < l_width; ++i)
std::cout << "_";
std::cout << std::endl;
}
///////////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::print_conv_footer()
{
static const unsigned int l_width = 102;
for (unsigned int i = 0; i < l_width; ++i)
std::cout << "_";
std::cout << std::endl;
std::cout << "Relative errors:" << std::endl
<< "Displacement:\t" << error_update.u / error_update_0.u << std::endl
<< "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl
<< "v / V_0:\t" << vol_current << " / " << vol_reference
<< std::endl;
}
/////////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::get_error_residual(Errors &error_residual)
{
BlockVector<double> error_res(dofs_per_block);
for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
if (!constraints.is_constrained(i))
error_res(i) = system_rhs(i);
error_residual.norm = error_res.l2_norm();
error_residual.u = error_res.block(u_dof).l2_norm();
}
////////////////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::get_error_update(const BlockVector<double> &newton_update,
Errors &error_update)
{
BlockVector<double> error_ud(dofs_per_block);
for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i)
if (!constraints.is_constrained(i))
error_ud(i) = newton_update(i);
error_update.norm = error_ud.l2_norm();
error_update.u = error_ud.block(u_dof).l2_norm();
}
// computed the nodal displacement vector
template <int dim>
BlockVector<double>
Solid<dim>::get_total_solution(const BlockVector<double> &solution_delta) const
{
BlockVector<double> solution_total(solution_n_iter);
solution_total += solution_delta;
return solution_total;
}
// next three functions compute elemental tangent stiffness and assemble the global tangent stiffness
template <int dim>
void Solid<dim>::assemble_system_tangent()
{
timer.enter_subsection("Assemble tangent matrix");
std::cout << " ASM_K " << std::flush;
tangent_matrix = 0.0;
const UpdateFlags uf_cell(update_gradients |
update_JxW_values);
PerTaskData_K per_task_data(dofs_per_cell);
ScratchData_K scratch_data(fe, qf_cell, uf_cell);
WorkStream::run(dof_handler_ref.begin_active(),
dof_handler_ref.end(),
*this,
&Solid::assemble_system_tangent_one_cell,
&Solid::copy_local_to_global_K,
scratch_data,
per_task_data);
timer.leave_subsection();
}
////////////////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::copy_local_to_global_K(const PerTaskData_K &data)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
tangent_matrix.add(data.local_dof_indices[i],
data.local_dof_indices[j],
data.cell_matrix(i, j));
}
////////////////////////////////////////////////////////////////////
template <int dim>
void
Solid<dim>::assemble_system_tangent_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_K &scratch,
PerTaskData_K &data)
{
data.reset();
scratch.reset();
scratch.fe_values_ref.reinit(cell);
cell->get_dof_indices(data.local_dof_indices);
PointHistory<dim> *lqph =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
const unsigned int k_group = fe.system_to_base_index(k).first.first;
if (k_group == u_dof)
{
scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point)
* F_inv;
scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]);
}
else
Assert(k_group <= u_dof, ExcInternalError());
}
}
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const Tensor<2, dim> tau = lqph[q_point].get_tau();
const SymmetricTensor<4, dim> Jc = lqph[q_point].get_Jc();
const std::vector<SymmetricTensor<2, dim> >
&symm_grad_Nx = scratch.symm_grad_Nx[q_point];
const std::vector<Tensor<2, dim> >
&grad_Nx = scratch.grad_Nx[q_point];
const double JxW = scratch.fe_values_ref.JxW(q_point);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const unsigned int component_i = fe.system_to_component_index(i).first;
const unsigned int i_group = fe.system_to_base_index(i).first.first;
for (unsigned int j = 0; j <= i; ++j)
{
const unsigned int component_j = fe.system_to_component_index(j).first;
const unsigned int j_group = fe.system_to_base_index(j).first.first;
if ((i_group == j_group) && (i_group == u_dof))
{
data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc
* symm_grad_Nx[j] * JxW;
if (component_i == component_j)
data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau
* grad_Nx[j][component_j] * JxW;
}
else
Assert((i_group <= u_dof) && (j_group <= u_dof),
ExcInternalError());
}
}
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
data.cell_matrix(i, j) = data.cell_matrix(j, i);
}
////////////////////////////////////////////////////
// next three functions calculate elemental residual and assemble them to obtain the global vector
template <int dim>
void Solid<dim>::assemble_system_rhs()
{
timer.enter_subsection("Assemble system right-hand side");
std::cout << " ASM_R " << std::flush;
system_rhs = 0.0;
const UpdateFlags uf_cell(update_values |
update_gradients |
update_JxW_values);
const UpdateFlags uf_face(update_values |
update_JxW_values);
PerTaskData_RHS per_task_data(dofs_per_cell);
ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face);
WorkStream::run(dof_handler_ref.begin_active(),
dof_handler_ref.end(),
*this,
&Solid::assemble_system_rhs_one_cell,
&Solid::copy_local_to_global_rhs,
scratch_data,
per_task_data);
timer.leave_subsection();
}
///////////////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::copy_local_to_global_rhs(const PerTaskData_RHS &data)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i);
}
///////////////////////////////////////////////////////////////////
template <int dim>
void
Solid<dim>::assemble_system_rhs_one_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_RHS &scratch,
PerTaskData_RHS &data)
{
data.reset();
scratch.reset();
scratch.fe_values_ref.reinit(cell);
cell->get_dof_indices(data.local_dof_indices);
PointHistory<dim> *lqph =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
const unsigned int k_group = fe.system_to_base_index(k).first.first;
if (k_group == u_dof)
scratch.symm_grad_Nx[q_point][k]
= symmetrize(scratch.fe_values_ref[u_fe].gradient(k, q_point)
* F_inv);
else
Assert(k_group <= u_dof, ExcInternalError());
}
}
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const SymmetricTensor<2, dim> tau = lqph[q_point].get_tau();
const std::vector<SymmetricTensor<2, dim> >
&symm_grad_Nx = scratch.symm_grad_Nx[q_point];
const double JxW = scratch.fe_values_ref.JxW(q_point);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const unsigned int i_group = fe.system_to_base_index(i).first.first;
if (i_group == u_dof)
data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
else
Assert(i_group <= u_dof, ExcInternalError());
}
}
////////////////
}
template <int dim>
void Solid<dim>::make_constraints(const int &it_nr)
{
std::cout << " CST " << std::flush;
if (it_nr > 1)
return;
constraints.clear();
const bool apply_dirichlet_bc = (it_nr == 0);
const FEValuesExtractors::Scalar x_displacement(0);
{
const int boundary_id = 1;
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
fe.component_mask(x_displacement));
else
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
fe.component_mask(x_displacement));
}
{ const double tol_boundary = 1.0e-12;
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler_ref.begin_active(),
endc = dof_handler_ref.end();
for (; cell!=endc; ++cell)
for (unsigned int v=0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
if ((std::abs(cell->vertex(v)[0] - 0.0) < tol_boundary) &&
(std::abs(cell->vertex(v)[1] -0.0)< tol_boundary) )
constraints.add_line(cell->vertex_dof_index(v, 1));
}
constraints.close();
}
////////////////////////////////////////////////////////////////
template <int dim>
std::pair<unsigned int, double>
Solid<dim>::solve_linear_system(BlockVector<double> &newton_update)
{
BlockVector<double> A(dofs_per_block);
BlockVector<double> B(dofs_per_block);
unsigned int lin_it = 0;
double lin_res = 0.0;
{
timer.enter_subsection("Linear solver");
std::cout << " SLV " << std::flush;
if (parameters.type_lin == "CG")
{
const int solver_its = tangent_matrix.block(u_dof, u_dof).m()
* parameters.max_iterations_lin;
const double tol_sol = parameters.tol_lin
* system_rhs.block(u_dof).l2_norm();
SolverControl solver_control(solver_its, tol_sol);
GrowingVectorMemory<Vector<double> > GVM;
SolverCG<Vector<double> > solver_CG(solver_control, GVM);
PreconditionSelector<SparseMatrix<double>, Vector<double> >
preconditioner (parameters.preconditioner_type,
parameters.preconditioner_relaxation);
preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
newton_update.block(u_dof),
system_rhs.block(u_dof),
preconditioner);
lin_it = solver_control.last_step();
lin_res = solver_control.last_value();
}
else if (parameters.type_lin == "Direct")
{
SparseDirectUMFPACK A_direct;
A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
A_direct.vmult(newton_update.block(u_dof), system_rhs.block(u_dof));
lin_it = 1;
lin_res = 0.0;
}
else
Assert (false, ExcMessage("Linear solver type not implemented"));
timer.leave_subsection();
}
constraints.distribute(newton_update);
return std::make_pair(lin_it, lin_res);
}
///////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::output_results() const
{
std::vector<std::string> solution_name(dim, "displacement");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(dim,
DataComponentInterpretation::component_is_part_of_vector);
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler_ref);
data_out.add_data_vector(solution_n,
solution_name,
DataOut<dim>::type_dof_data,
data_component_interpretation);
data_out.add_data_vector (dof_handler_eta, solution_etai, "etai");
Vector<double> soln(solution_n.size());
for (unsigned int i = 0; i < soln.size(); ++i)
soln(i) = solution_n(i);
MappingQEulerian<dim> q_mapping(degree, dof_handler_ref, soln);
data_out.build_patches(q_mapping, degree);
std::ostringstream filename;
filename << "solution-" << time.get_timestep() << ".vtk";
std::ofstream output(filename.str().c_str());
data_out.write_vtk(output);
}
}
/////////////////////////////////////////////////////////
int main ()
{
using namespace dealii;
using namespace PhaseField;
try
{
deallog.depth_console(0);
Solid<2> solid_3d("parameters.prm");
solid_3d.run();
}
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;
}
return 0;
}
parameters.prm
Description: Binary data
