Dear Jean,
Thanks very much for your help and support. Most of my problem is related
to applying axisymmetry boundary condition for a curved domain. I attached
my code for your consideration. According to the fact that this is my first
project with dealii, Please accept my lack of knowledge in this case. This
is a themo-elastic problem for a curved domain which for getting result, I
have fixed the straight line of my domain. However, I really do not know
why the Neumann boundary condition which is brought in the
assemble-system-rhs does not work properly, instead of putting force on the
curved domain, it seems that it is a tensile force. It should be noted that
I have changed the sign of my equation for this type of boundary condition,
however, it still has not worked properly. Your noteworthy comments on my
code can give me incredible experiences on working with dealii. Looking
forward to hearing from you.
Bests,
Benhour
On Monday, December 5, 2016 at 5:01:25 AM UTC-6, Jean-Paul Pelteret wrote:
>
> Dear Benhour,
>
>
>> These are 9 concepts, in 9 lines of text, all without any kind of
>> elaboration.
>> In other words, we have no real idea what you are doing, and equally
>> important, we have no idea what is actually wrong. When debugging
>> problems
>> like yours, you need to learn to reduce things to the minimal case that
>> shows
>> the problem, and to look for the *first sign* where something seems to be
>> going wrong.
>
>
> Prof Bangerth has succinctly stated why I'm having such difficulties in
> helping you.
>
> I, like him, would encourage you to take your code and strip out as many
> of the complexities as possible, and produce some test cases that convince
> you that each part of your code is working as expected. This is exactly how
> I build up any of my codes that tackle complex problems. For this you often
> need to test your implementation of each "feature" against a really simple
> problem that you can hand calculate a comparison result for, or otherwise
> some benchmark in the literature.
>
> Thats really the best advice we can give you in this instance. So, for
> example, coming back to your question of the implementation of the traction
> boundary condition: you could easily test this by considering a distributed
> load on a bar (inducing uni-axial tension) or perhaps the result shown in
> step-44 <https://www.dealii.org/8.4.1/doxygen/deal.II/step_44.html> or
> one of the more simple code-gallery examples
> <https://dealii.org/developer/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html>
> .
>
> Regards,
> Jean-Paul
>
--
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.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2000 - 2015 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 at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, University of Heidelberg, 2000
*/
// Just as in previous examples, we have to include several files of which the
// meaning has already been discussed:
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/base/std_cxx11/shared_ptr.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/compressed_sparsity_pattern.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/sparse_ilu.h>
#include <deal.II/lac/precondition_selector.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_dgp_monomial.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/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
namespace Melting
{
using namespace dealii;
namespace Parameters
{
/////////////////////////////////////////////////
struct MaterialVariables
{
double lambda0; // Bulk modulus of liquid phase
double lambda1; // Bulk modulus of solid phase
double mu1; // Shear modulus of solid phase (The shear modulus of liquid phase is assumed to be zero)
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void MaterialVariables::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Material properties");
{
prm.declare_entry("bulk modulus of liquid phase", "41.3e9", /* */
Patterns::Double(),
"bulk modulus of liquid phase");
prm.declare_entry("bulk modulus of solid phase", "71.1e9",
Patterns::Double(0.0),
"bulk modulus of solid phase");
prm.declare_entry("shear modulus of solid phase", "27.26e9",
Patterns::Double(),
"shear modulus of solid phase");
}
prm.leave_subsection();
}
void MaterialVariables::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Material properties");
{
lambda0 = prm.get_double("bulk modulus of liquid phase");
lambda1 = prm.get_double("bulk modulus of solid phase");
mu1 = prm.get_double("shear modulus of solid phase");
}
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", // changed from "1e-6", 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", // changed from "1.0e-9", 1.0e-15
Patterns::Double(0.0),
"Force residual tolerance");
prm.declare_entry("Tolerance displacement", "1.0e-12", // changed from "1.0e-6"
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 time_step; // It is equal to delta_t
double final_time; // It is equal to 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", "4.5e-10",
Patterns::Double(),
"End time");
prm.declare_entry("Time step size", "0.1e-12",
Patterns::Double(),
"Time step size");
}
prm.leave_subsection();
}
void Time::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Time");
{
final_time = prm.get_double("End time");
time_step = prm.get_double("Time step size");
}
prm.leave_subsection();
}
////////////////////////////////////////////////////////////
struct AllParameters :
public MaterialVariables,
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)
{
MaterialVariables::declare_parameters(prm);
LinearSolver::declare_parameters(prm);
NonlinearSolver::declare_parameters(prm);
Time::declare_parameters(prm);
}
void AllParameters::parse_parameters(ParameterHandler &prm)
{
MaterialVariables::parse_parameters(prm);
LinearSolver::parse_parameters(prm);
NonlinearSolver::parse_parameters(prm);
Time::parse_parameters(prm);
}
}
//////DEFINE SECOND ORDER IDENTITY TENSOR
template <int dim>
class StandardTensors
{
public:
static const SymmetricTensor<2, dim> I;
};
template <int dim>
const SymmetricTensor<2, dim>
StandardTensors<dim>::I = unit_symmetric_tensor<dim>();
//////////////////// DEFINE TIME STEP, CURRENT TIME ETC.
class Time
{
public:
Time (const double final_time,
const double time_step)
:
timestep(0),
time_current(0.0),
final_time(final_time),
time_step(time_step)
{}
virtual ~Time()
{}
double current() const
{
return time_current;
}
double end() const
{
return final_time;
}
double get_time_step() const
{
return time_step;
}
unsigned int get_timestep() const
{
return timestep;
}
void increment()
{
time_current += time_step;
++timestep;
}
private:
unsigned int timestep;
double time_current;
const double final_time;
const double time_step;
};
///////////////////////////////////////////////////////////////
template<int dim>
class Materials
{
public:
Materials (const double lambda0,
const double lambda1,
const double mu1
)
:
eV(0.0), // The trace of total strain
lambda00(lambda0),
lambda11(lambda1),
mu11(mu1)
{}
~Materials()
{}
void update_material_data (const SymmetricTensor<2, dim> &total_strain,
const double phi)
{
eV = trace(total_strain); /////// The total Volumetric Strain
mu = mu11 * phi;
lambda = (lambda00 + (lambda11 - lambda00) * phi) - 2.0/3.0 * mu;
}
SymmetricTensor<4, dim> get_Jc() const
{
return get_Jc_modulus();
}
protected:
double eV; // The total Volumetric Strain
double lambda00;
double lambda11;
double mu11;
double lambda;
double mu;
/////////////// compute the stress_strain_tensor
SymmetricTensor<4, dim> get_Jc_modulus () const
{
SymmetricTensor<4, dim> elasticity_stiffness_tensor;
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)
elasticity_stiffness_tensor[i][j][k][l] = (((i==k) && (j==l) ? mu : 0.0) +
((i==l) && (j==k) ? mu : 0.0) +
((i==j) && (k==l) ? lambda : 0.0));
return elasticity_stiffness_tensor;
}
};
///////////////////////// Update the quadrature point history
template<int dim>
class PointHistory
{
public:
PointHistory()
:
material(NULL),
elasStress(SymmetricTensor<2, dim>()),
Jc(SymmetricTensor<4, dim>())
{}
virtual ~PointHistory()
{
delete material;
material = NULL;
}
void setup_local_quadrature_point (const Parameters::AllParameters ¶meters)
{
material = new Materials<dim>(parameters.lambda0, parameters.lambda1,
parameters.mu1);
update_values(Tensor<2, dim>(), double ());
}
void update_values (const Tensor<2, dim> &Grad_u, const double eta)
{
// const SymmetricTensor<2, dim> total_strain = 0.5*(Grad_u + transpose(Grad_u));
const SymmetricTensor<2, dim> total_strain = symmetrize(Grad_u);
const double phi = std::pow(eta, 2) * (3 - 2 * eta);
dphi = 6 * eta * (1 - eta);
const Tensor<2, dim> T_transformation = (0.06*(1-phi)*Tensor<2, dim>(StandardTensors<dim>::I))/3.0;
const Tensor<2, dim> T_thermal = 640.52 * 3.032e-5 * Tensor<2, dim>(StandardTensors<dim>::I) -
3.67*(4.268e-5-(1.236e-5)*phi) * Tensor<2, dim>(StandardTensors<dim>::I);
Tensor<2, dim> tmp = T_transformation + T_thermal;
SymmetricTensor<2, dim> T_inelastic = symmetrize(tmp);
material->update_material_data(total_strain, phi);
T_elastic = total_strain - T_inelastic;
dev_T_elastic = T_elastic - (1/dim) * trace(T_elastic) * (Tensor<2, dim>(StandardTensors<dim>::I))/3.0;
Tensor<2, dim> tmpElas = (eV-T_inelastic*StandardTensors<dim>::I)*(41.3e9+(71.1e9-41.3e9)*phi)*
StandardTensors<dim>::I+2*27.3e9*phi*dev_T_elastic;
elasStress = symmetrize(tmpElas);
Jc = material->get_Jc(); // Extracting the elasticity modulus
get_eta_boundary_value = -2*((1.05-0.1292)+0.1292*phi) / 30e-9; // Extracting the boundary value for the mechanics
}
double get_eV() const
{
return material->get_eV();
}
const SymmetricTensor<2, dim> &get_elasStress() const
{
return elasStress;
}
const SymmetricTensor<4, dim> &get_Jc() const
{
return Jc;
}
const SymmetricTensor<2, dim> &elastic_strain() const
{
return T_elastic;
}
const SymmetricTensor<2, dim> &deviatoric_elastic_strain() const
{
return dev_T_elastic;
}
double get_eta_boundary_value_mechanic() const
{
return get_eta_boundary_value;
}
private:
Materials<dim> *material;
SymmetricTensor<2, dim> elasStress;
SymmetricTensor<4, dim> Jc;
SymmetricTensor<2, dim> T_elastic;
SymmetricTensor<2, dim> dev_T_elastic;
double get_eta_boundary_value;
double eV;
double dphi;
};
template <int dim>
class MeltingSolidification
{
public:
MeltingSolidification (const std::string &input_file);
virtual ~MeltingSolidification();
void run ();
private:
struct PerTaskData_K;
struct ScratchData_K;
struct PerTaskData_RHS;
struct ScratchData_RHS;
struct PerTaskData_UQPH;
struct ScratchData_UQPH;
void make_grid_and_dofs ();
void system_setup ();
//////////assemble the tangent stiffness for N-R problem for solving nodal displacements of Mechanical parts
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);
///////// Assemble system right hand side for the mechanics part or right hand side or residual
void assemble_system_right_hand_side ();
void assemble_system_right_hand_side_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_RHS &scratch,
PerTaskData_RHS &data);
void copy_local_to_global_right_hand_side(const PerTaskData_RHS &data);
/////// The tangent matrix for the right hand side of the G-L equation
void compute_nl_matrix (const Vector<double> &old_data,
const Vector<double> &new_data,
SparseMatrix<double> &nl_matrix) const;
/////// The Fh(Un,l) for the G-L equation, the residual or the right hand side for the G-L equation
void compute_nl_term (const Vector<double> &old_data,
const Vector<double> &new_eta,
Vector<double> &nl_term) const;
///////// Impose boundary condition for mechanics or displacement
void
make_constraints(const int &it_nr);
//////////This one is called before the first timestep to set up a pristine state for the history variables
void setup_quadrature_point_history ();
//////// Update the displacement for computing the stress
void update_quadrature_point_history_incremental (const BlockVector<double> &solution_delta);
void
update_quadrature_point_history_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*/)
{}
///////// solves the non-linear mechanics problem to compute displacements
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;
////// computes the values and the updated values of the order parameter
unsigned int
solve_eta ();
// void
// compute_eta_rhs (Vector<double> &force_eta) const;
void output_results () const;
Parameters::AllParameters parameters;
Triangulation<dim> triangulation;
Time time; // variable of type class 'Time'
TimerOutput timer;
std::vector<PointHistory<dim> > quadrature_point_history;
const unsigned int degree; // degree of polynomial of shape functions
const FESystem<dim> fe; // FE object for the displacement
DoFHandler<dim> dof_handler_ref; // DoF for Mechanic parts
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
// eta_dof = 1 //added
};
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 stiffness 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
// const unsigned int n_q_points_eta; // numbr of quadrature point for eta in the cell
// const unsigned int n_face_q_points_eta; // numbr of quadrature point for eta in the face
SparsityPattern sparsity_pattern_eta;
SparseMatrix<double> system_matrix; // system_matrix for G-L equation
SparseMatrix<double> mass_matrix; // mass matrix out of Ginxburg-Landau eqn
SparseMatrix<double> laplace_matrix; // Laplace matrix out of Ginxburg-Landau eqn
Vector<double> solution_delta_eta;
Vector<double> solution_eta; // solution vector for eta
Vector<double> old_solution_eta;
Vector<double> system_rhs_eta; // a vector used for solving eta (system_rhs in my previous code)
//
const double theta, DelK, laH, delgam,
landa, betta;
//////// 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();
// const unsigned int output_timestep_skip;
};
/////// defining the initial value for the order parameter
template <int dim>
class InitialValues : public Function<dim>
{
public:
InitialValues () : Function<dim>() {}
virtual double value(const Point<dim> &p,
const unsigned int /*component = 0*/) const;
};
/////////////////////////////////////////////////////
template <int dim>
double InitialValues<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
return tanh((sqrt((p[0]*p[0])+(p[1]*p[1]))-15e-9)/1e-9)+1;
// Assert (component == 0, ExcInternalError());
// return 0;
}
//////////////////////////////////////////////// class for boundary for the mechanics problem
template<int dim>
class BoundaryValuesM: public Function<dim>
{
public:
BoundaryValuesM() :
Function<dim>(dim + 1)
{}
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value(const Point<dim> &p, Vector<double> &value) const;
};
template<int dim>
double BoundaryValuesM<dim>::value(const Point<dim> &p,
const unsigned int component) const
{
Assert(component < this->n_components,
ExcIndexRange (component, 0, this->n_components));
Assert(p[0] >= -10, ExcLowerRange (p[0], 0)); //value of -10 is to permit some small numerical error moving nodes left of x=0; a << value is in fact sufficient
return -1e9;
}
template<int dim>
void BoundaryValuesM<dim>::vector_value(const Point<dim> &p,
Vector<double> &values) const
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = BoundaryValuesM<dim>::value(p, c);
}
///////////////////////////////////////////////////////////////////////////////////
template <int dim>
MeltingSolidification<dim>::MeltingSolidification (const std::string &input_file)
:
parameters(input_file),
triangulation(Triangulation<dim>::maximum_smoothing),
time(parameters.final_time, parameters.time_step),
timer(std::cout,
TimerOutput::summary,
TimerOutput::wall_times),
degree(2),
fe(FE_Q<dim>(2), dim), // For displacement
dof_handler_ref(triangulation),
dofs_per_cell (fe.dofs_per_cell),
u_fe(first_u_component),
dofs_per_block(n_blocks),
qf_cell(3),
qf_face(3),
n_q_points (qf_cell.size()),
n_q_points_f (qf_face.size()),
degree_eta(2),
fe_eta (2), // For the order parameter
dof_handler_eta (triangulation),
dofs_per_cell_eta(fe_eta.dofs_per_cell),
// n_q_points_eta(qf_cell.size()),
// n_face_q_points_eta(qf_cell.size()),
theta (0.5),
DelK (29.8e9),
laH(933.57e6),
delgam (0.1292),
landa (532.0),
betta (3.21e-10)
{
// determine_component_extractor();
}
/////////Destructor
template <int dim>
MeltingSolidification<dim>::~MeltingSolidification()
{
dof_handler_ref.clear();
dof_handler_eta.clear();
}
///////// The RUN
template <int dim>
void MeltingSolidification<dim>::run() // assemble system (step 25) comes in this section...
{
make_grid_and_dofs(); // generates the geometry and mesh
system_setup(); // sets up the system matrices
VectorTools::interpolate(dof_handler_eta, InitialValues<dim>(), solution_eta); //initial eta
output_results(); // prints output
/////// This is the zeroth iteration for compute the initial distorted solution
/////// Of the body due to arbitrary distribution of initial eta
BlockVector<double> solution_delta(dofs_per_block); // du (displacement vector)
solution_delta = 0.0;
solution_n_iter = solution_n; // transferring displacement soln to another vector which
//will be called from get_total_solution ()
solve_nonlinear_timestep(solution_delta);
solution_n += solution_delta; // obtaining the total displacement
output_results();
time.increment();
Vector<double> solution_delta_eta(solution_eta.size());
Vector<double> system_rhs_eta (solution_eta.size()); // the right hand side for the G-L equation
/////// computed actual time integration to update displacement and eta
while (time.current() <= time.end())
{
/////// first we solve for the order parameter
solution_delta_eta = 0.0;
///////First we assemble the Jacobian matrix F′h(Un,l), where Un,l is stored in the vector solution for convenience.
system_matrix.copy_from (mass_matrix);//system_matrix is the matrix that we want to invert it.
system_matrix.add (std::pow(parameters.time_step*theta*landa*betta,1), laplace_matrix);
SparseMatrix<double> tmp_matrix (sparsity_pattern_eta);
compute_nl_matrix (old_solution_eta, solution_eta, tmp_matrix);
system_matrix.add (-std::pow(parameters.time_step*theta*landa,1), tmp_matrix);
////Then, we compute the right-hand side vector −Fh(Un,l).
system_rhs_eta = 0;
tmp_matrix.copy_from (mass_matrix);
tmp_matrix.add (std::pow(parameters.time_step*theta*landa*betta,1), laplace_matrix);
Vector<double> tmp_vector (solution_eta.size());
tmp_matrix.vmult (tmp_vector, solution_eta);
system_rhs_eta += tmp_vector;
tmp_matrix.add(-std::pow(parameters.time_step*landa*betta, 1), laplace_matrix);
tmp_matrix.vmult (tmp_vector, old_solution_eta);
system_rhs_eta -= tmp_vector;
compute_nl_term (old_solution_eta, solution_eta, tmp_vector);
system_rhs_eta.add (std::pow(parameters.time_step*landa,1),tmp_vector);
system_rhs_eta *= -1;
////////////////////////////////////////// I add these lines for solving the order parameter by Newton-Raphson Method
std::cout << std::endl << "Timestep " << time.get_timestep() << " @ "
<< time.current() << "s" << std::endl;
unsigned int newton_iteration = 0;
for (; newton_iteration < parameters.max_iterations_NR;
++newton_iteration)
{
old_solution_eta = solution_eta;
double initial_rhs_norm = 0;
bool first_iteration = true;
do{
if (first_iteration == true)
initial_rhs_norm = system_rhs_eta.l2_norm();
const unsigned int newton_iteration = solve_eta();
solution_eta += solution_delta_eta; // Getting the total eta
if (first_iteration == true)
std::cout << " " << newton_iteration;
else
std::cout << '+' << newton_iteration;
first_iteration = false;
}
while (system_rhs_eta.l2_norm()>1e-6 * initial_rhs_norm);
std::cout << " CG iterations per nonlinear step."
<< std::endl;
////////////////////////////////////////////////////////////////////////////////////////////////////////////
// unsigned int timestep = 1;
// for (;time.current()<=time.end(); ++timestep)
// {
// old_solution_eta = solution_eta;
// std::cout << std::endl
// << "TimeStep #" << time.get_timestep() << "; "
// << "advancing to t= " << time.current() << "." << std::endl;
//
// double initial_rhs_norm = 0;
// bool first_iteration = true;
//
// do{
// if (first_iteration == true)
// initial_rhs_norm = system_rhs_eta.l2_norm();
// const unsigned int newton_iteration = solve_eta();
// solution_eta += solution_delta_eta; // Getting the total eta
// if (first_iteration == true)
// std::cout << " " << newton_iteration;
// else
// std::cout << '+' << newton_iteration;
// first_iteration = false;
// }
// while (system_rhs_eta.l2_norm()>1e-6 * initial_rhs_norm);
// std::cout << " CG iterations per nonlinear step."
// << std::endl;
output_results ();
// solve_eta(solution_delta_eta); // Solving for the order parameter
update_quadrature_point_history_incremental(solution_delta); // updating the quadrature point history
///////Solve for the updated displacements
solution_delta = 0.0;
solution_n_iter = solution_n;
solve_nonlinear_timestep(solution_delta);
std::cout << solution_delta.l2_norm () << std::endl;
solution_n += solution_delta;
output_results();
time.increment();
}
}
}
/////////Next three functions are used in computing and assembling the system tangent stiffness matrix
template <int dim>
struct MeltingSolidification<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 MeltingSolidification<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;
}
}
}
};
////////Next three functions are used in computing and assembling the system right hand side or the residual
template <int dim>
struct MeltingSolidification<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 MeltingSolidification<dim>::ScratchData_RHS
{
FEValues<dim> fe_values_ref;
FEFaceValues<dim> fe_face_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_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)),
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_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),
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;
}
}
}
};
/////////Next two functions are used in updating the quadrature point history with new stress, driving force,
//////// and other quantities
template <int dim>
struct MeltingSolidification<dim>::PerTaskData_UQPH
{
void reset()
{}
};
////////////////////////////////////////////////////
template <int dim>
struct MeltingSolidification<dim>::ScratchData_UQPH
{
const BlockVector<double> &solution_total;
std::vector<Tensor<2, dim> > solution_grads_u_total;
std::vector<double> solution_eta_total; //check if required
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_eta_total (qf_cell.size()),
fe_values_ref(fe_cell, qf_cell, uf_cell),
fe_values_eta(fe_cell_eta, qf_cell, uf_cell_eta) // check if qf_cell_eta is required
{}
ScratchData_UQPH(const ScratchData_UQPH &rhs)
:
solution_total(rhs.solution_total),
solution_grads_u_total(rhs.solution_grads_u_total),
solution_eta_total(rhs.solution_eta_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_eta_total[q] = 0.0;
}
}
};
///////// The tangent matrix (nonlinear matrix) for computing the order parameter (equal to assemble system tangent)
template <int dim>
void MeltingSolidification<dim>::compute_nl_matrix (const Vector<double> &old_data,
const Vector<double> &new_data,
SparseMatrix<double> &nl_matrix) const
{
QGauss<dim> quadrature_formula_eta (3);
FEValues<dim> fe_values_eta (fe_eta, quadrature_formula_eta,
update_values | update_JxW_values | update_quadrature_points);
const unsigned int dofs_per_cell_eta = fe_eta.dofs_per_cell;
const unsigned int n_q_points_eta = quadrature_formula_eta.size();
FullMatrix<double> local_nl_matrix (dofs_per_cell_eta, dofs_per_cell_eta);
std::vector<types::global_dof_index> local_dof_indices_eta (dofs_per_cell_eta);
std::vector<double> old_data_values (n_q_points_eta);
std::vector<double> new_data_values (n_q_points_eta);
typename DoFHandler<dim>::active_cell_iterator
cell_eta = dof_handler_eta.begin_active(),
endc = dof_handler_eta.end();
for (; cell_eta!=endc; ++cell_eta)
{
local_nl_matrix = 0;
fe_values_eta.reinit (cell_eta);
fe_values_eta.get_function_values (old_data, old_data_values);
fe_values_eta.get_function_values (new_data, new_data_values);
PointHistory<dim> *local_quadrature_point_history =
reinterpret_cast<PointHistory<dim>*>(cell_eta->user_pointer());
for (unsigned int q_point=0; q_point<n_q_points_eta; ++q_point)
{
const SymmetricTensor<2, dim> elastic_stress = local_quadrature_point_history[q_point].get_elasStress();
const SymmetricTensor<2, dim> elastic_strain = local_quadrature_point_history[q_point].elastic_strain();
const SymmetricTensor<2, dim> deviatoric_strain = local_quadrature_point_history[q_point].deviatoric_elastic_strain();
for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
for (unsigned int j=0; j<dofs_per_cell_eta; ++j)
/////////////////////////////////////////////////////////////////////////////////////
////////////// I add a new term for nonlinear matrix according to my problem
local_nl_matrix(i,j) += ((12*560.13e6*(std::pow(theta * new_data_values[q_point] +
(1-theta) * old_data_values[q_point],2)-12*560.13e6*std::pow(theta *
new_data_values[q_point] +(1-theta) * old_data_values[q_point],1)+2*560.13e6+
6*(0.5*29.8e9
*(elastic_strain*elastic_strain)+27.26e9
*(deviatoric_strain*deviatoric_strain)
-3.6696e6+0.02*elastic_stress*StandardTensors<dim>::I / 3.0-
elastic_stress*StandardTensors<dim>::I*1.236e-5*3.6696)-12*
(0.5*29.8e9
*(elastic_strain*elastic_strain)
+27.26e9
*(deviatoric_strain*deviatoric_strain)
-
3.6696e6+0.02*elastic_stress*StandardTensors<dim>::I / 3.0 -
elastic_stress*StandardTensors<dim>::I*1.236e-5*3.6696)*
std::pow(theta*new_data_values[q_point]+(1-theta)*old_data_values[q_point],1)))*
fe_values_eta.shape_value (i, q_point) *
fe_values_eta.shape_value (j, q_point) *
fe_values_eta.JxW (q_point));
/////////////////////////////////////////////////////////////////////////////////////////////////////////
/// This is equal to copy local to global tangent matrix
cell_eta->get_dof_indices (local_dof_indices_eta);
for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
for (unsigned int j=0; j<dofs_per_cell_eta; ++j)
nl_matrix.add(local_dof_indices_eta[i], local_dof_indices_eta[j],
local_nl_matrix(i,j));
}
}
}
///////////Computes the right hand side (non-linear term) in the Ginzburg-Landau equation (equal to assemble system rhs)
template <int dim>
void MeltingSolidification<dim>::compute_nl_term (const Vector<double> &old_data,
const Vector<double> &new_data,
Vector<double> &nl_term) const
{
nl_term = 0;
const QGauss<dim> quadrature_formula_eta (3);
FEValues<dim> fe_values_eta (fe_eta, quadrature_formula_eta,
update_values |
update_JxW_values |
update_quadrature_points);
QGauss<dim-1> face_quadrature_formula_eta(3);
FEFaceValues<dim> fe_face_values_eta (fe_eta, face_quadrature_formula_eta,
update_values | update_gradients | update_JxW_values);
const unsigned int dofs_per_cell_eta = fe_eta.dofs_per_cell;
const unsigned int n_face_q_points_eta = face_quadrature_formula_eta.size();
const unsigned int n_q_points_eta = quadrature_formula_eta.size();
Vector<double> local_nl_term (dofs_per_cell_eta);
std::vector<types::global_dof_index> local_dof_indices_eta (dofs_per_cell_eta);
std::vector<double> old_data_values (n_q_points_eta);
std::vector<double> new_data_values (n_q_points_eta);
typename DoFHandler<dim>::active_cell_iterator
cell_eta = dof_handler_eta.begin_active(),
endc = dof_handler_eta.end();
for (; cell_eta!=endc; ++cell_eta)
{
local_nl_term = 0;
fe_values_eta.reinit (cell_eta);
fe_values_eta.get_function_values (old_data, old_data_values);
fe_values_eta.get_function_values (new_data, new_data_values);
PointHistory<dim> *local_quadrature_point_history =
reinterpret_cast<PointHistory<dim>*>(cell_eta->user_pointer());
for (unsigned int q_point=0; q_point<n_q_points_eta; ++q_point)
{
const SymmetricTensor<2, dim> elastic_stress = local_quadrature_point_history[q_point].get_elasStress();
const SymmetricTensor<2, dim> elastic_strain = local_quadrature_point_history[q_point].elastic_strain();
const SymmetricTensor<2, dim> deviatoric_strain = local_quadrature_point_history[q_point].deviatoric_elastic_strain();
for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
{
/////////////////////////////////////////////////////////////////////
// I add new right hand side for the nl_term
local_nl_term(i) += ((4*560.13e6*std::pow(theta * new_data_values[q_point] +
(1-theta) * old_data_values[q_point],3)-6*560.13e6*
std::pow(theta * new_data_values[q_point] +
(1-theta) * old_data_values[q_point],2)+2*560.13e6*
std::pow(theta * new_data_values[q_point] + (1-theta) * old_data_values[q_point],1)+
6*(0.5*29.8e9
*(elastic_strain * elastic_strain)
+27.26e9
*(deviatoric_strain*deviatoric_strain)
-3.6696e6+0.02*elastic_stress*StandardTensors<dim>::I / 3.0 -
elastic_stress*StandardTensors<dim>::I*
1.236e-5*3.6696)*std::pow(theta*new_data_values[q_point]+(1-theta)*old_data_values[q_point],1)-
6*(0.5*29.8e9
*(elastic_strain*elastic_strain)
+27.26e9
*(deviatoric_strain*deviatoric_strain)
-3.6696e6+0.02*elastic_stress*StandardTensors<dim>::I / 3.0 -
elastic_stress*StandardTensors<dim>::I*
1.236e-5*3.6696)*std::pow(theta*new_data_values[q_point]+(1-theta)*old_data_values[q_point],2))*
fe_values_eta.shape_value (i, q_point) *
fe_values_eta.JxW (q_point));
}
}
////////////// I add these lines for computing the nonlinear term due to boundary conditions
for (unsigned int face_number=0; face_number<GeometryInfo<dim>::faces_per_cell; ++face_number)
if (cell_eta->face(face_number)->at_boundary()
&&
(cell_eta->face(face_number)->boundary_id() == 0))
{
fe_face_values_eta.reinit (cell_eta, face_number);
for (unsigned int f_q_point_eta=0; f_q_point_eta<n_face_q_points_eta; ++f_q_point_eta)
{
for (unsigned int i = 0; i<dofs_per_cell_eta; ++i)
local_nl_term(i) -= ((6*theta*delgam*landa*std::pow(theta * new_data_values[f_q_point_eta] +
(1-theta) * old_data_values[f_q_point_eta],2) -
6*delgam*landa*theta*std::pow(theta * new_data_values[f_q_point_eta]
+(1-theta) * old_data_values[f_q_point_eta],1))*
fe_face_values_eta.shape_value(i, f_q_point_eta) *
fe_face_values_eta.JxW(f_q_point_eta));
}
}
/////////////////////////////////////////////////////////////////////////////////////////////////
cell_eta->get_dof_indices (local_dof_indices_eta);
for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
nl_term(local_dof_indices_eta[i]) += local_nl_term(i);
}
}
//////////// Assemble system for the order parameter and G-L equation
//template<int dim>
//void MeltingSolidification<dim>::assemble_system_right_hand_side_eta ()
//{
// QGauss<dim> quadrature_formula_eta(3); //(fe.degree+1);
////////// I add this for computing the boundary condition on the face of the outer curved edge
// QGauss<dim-1> face_quadrature_formula_eta(3);
//
// const unsigned int dofs_per_cell_eta = fe_eta.dofs_per_cell;
// const unsigned int n_q_points_eta = quadrature_formula_eta.size();
//
///////// I add this for computing the number of quadrature points in the face
// const unsigned int n_face_q_points_eta = face_quadrature_formula_eta.size();
// FullMatrix<double> cell_mass_matrix(dofs_per_cell, dofs_per_cell);
// FullMatrix<double> cell_laplace_matrix(dofs_per_cell, dofs_per_cell);
// Vector<double> cell_rhs_eta (dofs_per_cell_eta);
//
// std::vector<types::global_dof_index> local_dof_indices_eta(dofs_per_cell_eta);
//
// FEValues<dim> fe_values_eta(fe_eta, quadrature_formula_eta,
// update_values | update_gradients | update_JxW_values);
//
///////// I add this for computing the boundary condition on the outer surface of the surface
// FEFaceValues<dim> fe_face_values_eta (fe_eta, face_quadrature_formula_eta,
// update_values | update_JxW_values);
// const Vector<double> old_data;
// const Vector<double> new_data;
// std::vector<double> old_data_values (n_q_points_eta);
// std::vector<double> new_data_values (n_q_points_eta);
//
// typename DoFHandler<dim>::active_cell_iterator
// cell_eta = dof_handler_eta.begin_active(),
// endc = dof_handler_eta.end();
// for (; cell_eta!=endc; ++cell_eta)
// {
// fe_values_eta.reinit(cell_eta);
// cell_mass_matrix = 0;
// cell_laplace_matrix = 0;
// cell_rhs_eta = 0;
// fe_values_eta.get_function_values (old_data, old_data_values);
// fe_values_eta.get_function_values (new_data, new_data_values);
//
// for (unsigned int q_index=0; q_index<n_q_points_eta; ++q_index)
// {
// for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
// for (unsigned int j=0; j<dofs_per_cell_eta; ++j)
// cell_mass_matrix(i,j) += (fe_values_eta.shape_value(i, q_index) *
// fe_values_eta.shape_value(j, q_index) *
// fe_values_eta.JxW(q_index));
//
// for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
// for (unsigned int j=0; j<dofs_per_cell_eta; ++j)
// cell_laplace_matrix(i,j) += (fe_values_eta.shape_grad(i, q_index) *
// fe_values_eta.shape_grad(j, q_index) *
// fe_values_eta.JxW(q_index));
//
//
// for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
// cell_rhs_eta(i) += (fe_values_eta.shape_value (i, q_index) *
// 1 *
// fe_values_eta.JxW (q_index));
// }
//
///////// I add these lines for computing the nonlinear term due to boundary conditions
//
// for (unsigned int face_number=0; face_number<GeometryInfo<dim>::faces_per_cell; ++face_number)
// if (cell_eta->face(face_number)->at_boundary()
// &&
// cell_eta->face(face_number)->boundary_id() == 2)
// {
// fe_face_values_eta.reinit (cell_eta, face_number);
// for (unsigned int q_point=0; q_point<n_face_q_points_eta; ++q_point)
// {
// for (unsigned int i = 0; i<dofs_per_cell_eta; ++i)
// cell_rhs_eta(i) += ((6*theta*delgam*landa*std::pow(theta * new_data_values[n_face_q_points_eta] +
// (1-theta) * old_data_values[n_face_q_points_eta],2) -
// 6*delgam*landa*theta*std::pow(theta * new_data_values[n_face_q_points_eta]
// +(1-theta) * old_data_values[n_face_q_points_eta],1))*
// fe_face_values_eta.shape_value(i, n_face_q_points_eta) *
// fe_face_values_eta.JxW(n_face_q_points_eta));
// }
// }
//
// cell_eta->get_dof_indices(local_dof_indices_eta);
// for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
// for (unsigned int j=0; j<dofs_per_cell_eta; ++j)
// mass_matrix.add(local_dof_indices_eta[i],
// local_dof_indices_eta[j],
// cell_mass_matrix(i,j));
// for (unsigned int i=0; i<dofs_per_cell; ++i)
// for (unsigned int j=0; j<dofs_per_cell; ++j)
// laplace_matrix.add(local_dof_indices_eta[i],
// local_dof_indices_eta[j],
// cell_laplace_matrix(i,j));
//
// for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
// system_rhs_eta(local_dof_indices_eta[i]) += cell_rhs_eta(i);
// }
//
//}
template <int dim>
void MeltingSolidification<dim>::make_grid_and_dofs ()
{
////////////////////////// I add this section for changing my boundary condition and geometry
// const Point<2> center;
// const double inner_radius = 15e-9,
// outer_radius = 30e-9;
// GridGenerator::quarter_hyper_shell (triangulation,
// center, inner_radius, outer_radius,
// 10, true);
// static const SphericalManifold<2> boundary_description(center);
// triangulation.set_manifold (0, boundary_description);
// Triangulation<2>::active_cell_iterator
// cell = triangulation.begin_active(),
// endc = triangulation.end();
// for (; cell!=endc; ++cell)
// cell->set_all_manifold_ids (0);
// triangulation.refine_global (4);
//////////////////////////////////////////////////////////////////////////////
Point <dim> center;
const double radius = 15e-9;
GridGenerator::half_hyper_ball(triangulation,
center,
radius);
static const SphericalManifold<dim> manifold_description(center);
triangulation.set_all_manifold_ids_on_boundary(0);
triangulation.set_manifold (0, manifold_description);
//
//
const Point<dim> mesh_center;
for (typename Triangulation<dim>::active_cell_iterator cell = triangulation.begin_active();
cell != triangulation.end(); ++cell)
if (mesh_center.distance (cell->center()) > cell->diameter()/10e9)
cell->set_all_manifold_ids (0);
GridTools::copy_boundary_to_manifold_id(triangulation);
triangulation.refine_global (5);
// const double angle = numbers::PI / 2.0;
// GridTools::rotate (angle, triangulation);
// const double tol_boundary = 1e-3;
// for (typename Triangulation<dim> ::active_cell_iterator
// cell=triangulation.begin_active();
// cell!=triangulation.end(); ++cell) // loop over all cells
// {
// for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) // loop over all vertices
// {
// if (cell->face(face)->at_boundary())
// {
// const Point<dim> face_center = cell->face(face)->center();
// if (std::abs(std::sqrt(face_center[0]*face_center[0]+face_center[1]*face_center[1])
// - radius) < tol_boundary)
// cell->face(face)->set_boundary_id (2); // faces on the outer curved edge of the domain...
// else if (std::abs(face_center[0]-0.0) < tol_boundary) // -x-axis
// cell->face(face)->set_boundary_id (1);
// else if (std::abs(face_center[0] - 30e-9) < tol_boundary) // +x-axis
// cell->face(face)->set_boundary_id(0);
//// else if (std::abs(face_center[1] - radius) < tol_boundary)
//// cell->face(face)->set_boundary_id (1); // faces on the r-axis
//// else if (std::abs(face_center[0] - radius)< tol_boundary)
//// cell->face(face)->set_boundary_id (0); // faces on the z-axis
// }
// }
// }
}
////////// Implying the setup system
template <int dim>
void MeltingSolidification<dim>::system_setup() // Implying the degrees of freedom
{
timer.enter_subsection("Setup system");
std::vector<unsigned int> block_component(n_components, u_dof); // Displacement
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();
{
system_matrix.clear ();
mass_matrix.clear ();
laplace_matrix.clear ();
DynamicSparsityPattern dsp (n_dofs_eta, n_dofs_eta);
DoFTools::make_sparsity_pattern (dof_handler_eta, dsp);
sparsity_pattern_eta.copy_from(dsp);
system_matrix.reinit (sparsity_pattern_eta);
mass_matrix.reinit (sparsity_pattern_eta);
laplace_matrix.reinit (sparsity_pattern_eta);
////////Creating the mass matrix
MatrixCreator::create_mass_matrix(dof_handler_eta,
QGauss<dim>(3),
mass_matrix);
///////Creating the Laplace matrix
MatrixCreator::create_laplace_matrix(dof_handler_eta,
QGauss<dim>(3),
laplace_matrix);
solution_eta.reinit(dof_handler_eta.n_dofs());
solution_delta_eta.reinit (dof_handler_eta.n_dofs());
old_solution_eta.reinit (dof_handler_eta.n_dofs());
system_rhs_eta.reinit(dof_handler_eta.n_dofs());; // system_rhs for the step25 or temp vector
}
// Setting up the quadrature point history
setup_quadrature_point_history();
timer.leave_subsection();
}
///////// Setup quadrature point history
template <int dim>
void MeltingSolidification<dim>::setup_quadrature_point_history()
{
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> *local_quadrature_point_history =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
Assert(local_quadrature_point_history >= &quadrature_point_history.front(), ExcInternalError());
Assert(local_quadrature_point_history <= &quadrature_point_history.back(), ExcInternalError());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
local_quadrature_point_history[q_point].setup_local_quadrature_point(parameters);
}
}
///////Updates the data at all quadrature points over a loop run by WorkStream::run
template <int dim>
void MeltingSolidification<dim>::update_quadrature_point_history_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);
PerTaskData_UQPH per_task_data_UQPH;
ScratchData_UQPH scratch_data_UQPH(fe, uf_UQPH, fe_eta, uf_UQPH_eta, qf_cell, solution_total); // I should add fe_eta
// and uf_UQPH_eta
WorkStream::run(dof_handler_ref.begin_active(),
dof_handler_ref.end(),
*this,
&MeltingSolidification::update_quadrature_point_history_incremental_one_cell,
&MeltingSolidification::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
MeltingSolidification<dim>::update_quadrature_point_history_incremental_one_cell(const typename
DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_UQPH &scratch,
PerTaskData_UQPH &/*data*/)
{
PointHistory<dim> *local_quadrature_point_history =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
Assert(local_quadrature_point_history >= &quadrature_point_history.front(), ExcInternalError());
Assert(local_quadrature_point_history <= &quadrature_point_history.back(), ExcInternalError());
Assert(scratch.solution_grads_u_total.size() == n_q_points,
ExcInternalError());
Assert(scratch.solution_eta_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); //check fe_values_eta to
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_eta, scratch.solution_eta_total// chech if scratch is needed
);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
local_quadrature_point_history[q_point].update_values(scratch.solution_grads_u_total[q_point],
scratch.solution_eta_total[q_point]);
}
/////Solving for the order parameter
template <int dim>
unsigned int
MeltingSolidification<dim>::solve_eta ()
{
SolverControl solver_control (1000, 1e-12*system_rhs.l2_norm());
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve (system_matrix, solution_delta_eta, system_rhs_eta,
preconditioner);
return solver_control.last_step();
}
/////// Solving the N-R problem for displacement
template <int dim>
void
MeltingSolidification<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_right_hand_side(); // 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_quadrature_point_history_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 MeltingSolidification<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 MeltingSolidification<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;
}
//////////////////
template <int dim>
void MeltingSolidification<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 MeltingSolidification<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>
MeltingSolidification<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 MeltingSolidification<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,
&MeltingSolidification::assemble_system_tangent_one_cell,
&MeltingSolidification::copy_local_to_global_K,
scratch_data,
per_task_data);
timer.leave_subsection();
}
///////////////////
template <int dim>
void MeltingSolidification<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));
}
/////////////////////local cell stiffness matrix
template <int dim>
void
MeltingSolidification<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> *local_quadrature_point_history =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
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);
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 SymmetricTensor<4, dim> Jc = local_quadrature_point_history[q_point].get_Jc();
const double JxW = scratch.fe_values_ref.JxW(q_point);
const std::vector<SymmetricTensor<2, dim> > &symm_grad_Nx = scratch.symm_grad_Nx[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;
for (unsigned int j = 0; j <= i; ++j)
{
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 // The material contribution:
* symm_grad_Nx[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 MeltingSolidification<dim>::assemble_system_right_hand_side()
{
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 |
update_normal_vectors);
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,
&MeltingSolidification::assemble_system_right_hand_side_one_cell,
&MeltingSolidification::copy_local_to_global_right_hand_side,
scratch_data,
per_task_data);
timer.leave_subsection();
}
///////////////////////////////////////////////////////////////////
template <int dim>
void MeltingSolidification<dim>::copy_local_to_global_right_hand_side(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
MeltingSolidification<dim>::assemble_system_right_hand_side_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> *local_quadrature_point_history =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
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));
else
Assert(k_group <= u_dof, ExcInternalError());
}
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
const SymmetricTensor<2, dim> elasStress = local_quadrature_point_history[q_point].get_elasStress();
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] * elasStress) * JxW;
else
Assert(i_group <= u_dof, ExcInternalError());
}
}
// Implying the mechanical boundary value to the outer surface of the domain
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary() == true
&& cell->face(face)->boundary_id() == 0)
{
scratch.fe_face_values_ref.reinit(cell, face);
for (unsigned int f_q_point=0; f_q_point<n_q_points_f; ++f_q_point)
{
const double time_ramp = (time.current() / time.end());
const double eta_orderM = local_quadrature_point_history[f_q_point].get_eta_boundary_value_mechanic() * time_ramp;
static const Tensor<1, dim> normals = scratch.fe_face_values_ref.normal_vector(f_q_point);
const Tensor<1, dim> neumann_value = eta_orderM * normals;
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)
{
const unsigned int component_i = fe.system_to_component_index(i).first;
const double JxW = scratch.fe_face_values_ref.JxW(f_q_point);
data.cell_rhs(i) += neumann_value[component_i] * scratch.fe_face_values_ref.shape_value(i, f_q_point) * JxW;
}
}
}
}
}
///////// Imposing the constraints for the mechanics problem, applying no_normal_flux_boundaries for boundary id 1
template <int dim>
void MeltingSolidification<dim>::make_constraints(const int &it_nr)
{
std::cout << " CST " << std::flush;
if (it_nr > 1)
return;
constraints.clear();
const bool apply_axisymmetry_bc = (it_nr == 0);
const FEValuesExtractors::Scalar x_displacement(0);
const FEValuesExtractors::Scalar y_displacement(1);
const FEValuesExtractors::Scalar z_displacement(2);
{
const int boundary_id = 1;
if (apply_axisymmetry_bc == true)
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
fe.component_mask(x_displacement) |
fe.component_mask(y_displacement)
);
else
VectorTools::interpolate_boundary_values(dof_handler_ref,
boundary_id,
ZeroFunction<dim>(n_components),
constraints,
fe.component_mask(x_displacement) |
fe.component_mask(y_displacement)
);
}
// {
// const int boundary_id = 1;
// if (apply_axisymmetry_bc == true)
// VectorTools::interpolate_boundary_values(dof_handler_ref,
// boundary_id,
// ZeroFunction<dim>(n_components),
// constraints,
// fe.component_mask(y_displacement)
// );
// else
// VectorTools::interpolate_boundary_values(dof_handler_ref,
// boundary_id,
// ZeroFunction<dim>(n_components),
// constraints,
// fe.component_mask(y_displacement)
// );
// }
///////////////////////////////////////////////////////////////////////////////
// {
// std::cout << " CST " << std::flush;
// if (it_nr > 1)
// return;
// constraints.clear();
// std::vector<bool> component_maskP(dim + 1, false);
// component_maskP[dim] = true;
// DoFTools::make_hanging_node_constraints(dof_handler_ref, constraints);
// VectorTools::interpolate_boundary_values(dof_handler_ref, 2,
// BoundaryValuesM<dim>(), constraints, component_maskP);
// }
// {
// std::cout << " CST " << std::flush;
// if (it_nr > 1)
// return;
// std::set<unsigned char> no_normal_flux_boundaries;
// no_normal_flux_boundaries.insert(1);
// VectorTools::compute_no_normal_flux_constraints(dof_handler_ref, 0,
// no_normal_flux_boundaries, constraints);
// }
// std::cout << " CST " << std::flush;
// if (it_nr > 1)
// return;
// constraints.clear();
// const bool apply_axialsymmetry_bc = (it_nr == 0);
// const FEValuesExtractors::Scalar r_displacement(0);
// const FEValuesExtractors::Scalar phi_displacement(1);
// const FEValuesExtractors::Scalar z_displacement(2);
//
// { // Implying n.u=0 for the r-axis of my problem
// std::set<unsigned char> no_normal_flux_boundaries;
// no_normal_flux_boundaries.insert (1);
// if (apply_axialsymmetry_bc == true)
// VectorTools::compute_no_normal_flux_constraints(dof_handler_ref,
// 0,
// no_normal_flux_boundaries,
// constraints
// );
// else
// VectorTools::compute_no_normal_flux_constraints(dof_handler_ref,
// 0,
// no_normal_flux_boundaries,
// constraints
// );
// }
// {
// if (apply_axisymmetry_bc == true)
// {
// typename DoFHandler<dim>::active_cell_iterator
// cell = dof_handler_ref.begin_active(),
// endc = dof_handler_ref.end();
// for (;cell!=endc; ++cell)
// {
// for (unsigned int vertex_index = 0; vertex_index < GeometryInfo<dim>::vertices_per_cell; ++vertex_index)
// {
// if ((cell->vertex(vertex_index)(0)==0.0) && (cell->vertex(vertex_index)(1)==0)
// // (cell->vertex(vertex_index)(2) == 0.5e-9)
// )
// {
// constraints.add_line(cell->vertex_dof_index(vertex_index, 0));
// constraints.add_line(cell->vertex_dof_index(vertex_index, 1));
// }
// }
// }
// }
// }
constraints.close();
}
////////////////////////////
template <int dim>
std::pair<unsigned int, double>
MeltingSolidification<dim>::solve_linear_system(BlockVector<double> &newton_update)
{
timer.enter_subsection("Linear solver");
BlockVector<double> A(dofs_per_block);
BlockVector<double> B(dofs_per_block);
unsigned int lin_it = 0;
double lin_res = 0.0;
{
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 MeltingSolidification<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_eta, "eta");
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_time_step() << ".vtk";
std::ofstream output(filename.str().c_str());
data_out.write_vtk(output);
}
}
int main (int argc, char *argv[])
{
using namespace dealii;
using namespace Melting;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, dealii::numbers::invalid_unsigned_int);
try
{
deallog.depth_console(0);
MeltingSolidification<2> melting_solidification_2d ("parameters.prm");
melting_solidification_2d.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;
}
