Dear All,
I want to solve a thermoelastic problem in cylindrical coordinate system.
My domain is a circle in 2D or Sphere in 3D. I am following step 44 with
considering the fact that I am solving for small strain. Because I have
curved domain and I should Impose a boundary condition that change with
temperature on the curved side of my domain, I understood that I should use
MappingQ because of the fact that my whole domain includes curved edge. For
this case, I have some problem with defining elasticity matrices such as
thermal, initial and elastic strain along with stresses. I have attached my
code for your consideration. It would be very kind of you if you take a
glance at it and tell me what is wrong defining my matrices.
On Tuesday, October 18, 2016 at 2:27:17 PM UTC-5, Jean-Paul Pelteret wrote:
>
> Dear Benhour,
>
> I'd like to help you but I really don't understand your question. Since
> its templated on the [space] dimension, the PointHistory class that you've
> displayed here is potentially suitable for either codimension 0 problems
> (dependent on your implementation of a constitutive law, of course, and
> assuming that you swap all instances of "dim" for the template parameter
> "spacedim"). Can you please try to explain a little more clearly what you
> are trying to do and what the issue is, and give some more context to the
> code you've quoted?
>
> - Are you solving your problem in the "standard" setting (either dim=1,2
> or 3), or the codimension 1 case (e.g. a 2-d manifold embedded in 3d
> space)? Or does one of your field variables only correspond to the
> codimension 1 case?
> - Are you referring to some structures that store, for example, the
> kinematic variables for both fields as well as the tangents resulting from
> linearisation of the coupled problem?
> - Are you solving your problem using a staggered scheme? I ask because you
> store neither any of the kinematic data nor tangents related to the thermal
> problem or coupling.
>
> Have you examined the many tutorials that solve problem with multiple
> fields/components? They may not be multi-physics *per se*, but will give
> you some insights on how to set up your problem.
>
> Regards,
> J-P
>
> On Tuesday, October 18, 2016 at 6:31:30 PM UTC+2, benhour....@gmail.com
> <javascript:> wrote:
>>
>> Dear All,
>> I am solving an elastic problem coupled with heat equation transfer
>> (thermal-elasticity). I should solve the problem in a curved domain (circle
>> in 2d and sphere in 3d) . I am developing step 44 in this case. I have to
>> define some symmetric tensors. I have attached some parts of my code that i
>> have problem with. I have two unknown variables. one for displacement
>> (vector) and another for other parameter(scalar). I do really appreciate
>> your kindness if you let me know how I can define symmetric tensor or in
>> general tensor for this case. In addition, should I define DofHandler<dim,
>> spacedim> for both scalar and vector valued unknown? My main problem is how
>> to define tensor in different dim and spacedim problem.
>>
>> template <int spacedim>
>> class PointHistory
>> {
>> public:
>> PointHistory()
>> :
>> material(NULL),
>> F_inv(StandardTensors<spacedim>::I),
>> elasStress(SymmetricTensor<2, spacedim>()),
>> totalStress(SymmetricTensor<2, spacedim>()),
>> Jc(SymmetricTensor<4, spacedim>())
>> {}
>>
>> virtual ~PointHistory()
>> {
>> delete material;
>> material = NULL;
>> }
>>
>> void setup_lqp (const Parameters::AllParameters ¶meters)
>> {
>> material = new
>> Material_Variables<spacedim>(parameters.lambda0,
>> parameters.mu0, parameters.lambda1, parameters.mu1,
>> parameters.chi, parameters.H,
>> parameters.thetain, parameters.thetacr, parameters.thetaeq,
>> parameters.theta, parameters.gammasv);
>>
>> update_values(Tensor<2, dim>(), double ());
>> }
>>
>> void update_values (const Tensor<2, dim> &Grad_u_n, const double
>> eta
>> )
>> {
>> const Tensor<2, dim> eps = 0.5 * (Grad_u_n +
>> transpose(Grad_u_n)); // Total
>> Strain
>>
>>
>> const double phi = std::pow(eta, 2) * (3 - 2 * eta);
>> // Function
>> of order parameter
>> dphi = 6 * eta * (1 - eta);
>>
>> // Derivative of the interpolation function
>> const double doubleWellBarrier = 2*eta + 4 * std::pow(eta, 3)
>> - 6*std::pow(eta, 2);
>> const double etaBVG = -6 * chi1 * delg1 * (1+eps_V) * dphi;
>> // Boundary value for
>> G-L equation
>> const double etaBVM = 2 * gam1 / 30e-9; // Boundary value for
>> Mechanics
>>
>> // Thermal Strain Tensor
>> const SymmetricTensor<2, dim> epsThermal = alphas * (thetaeq1 -
>> thetain1) * StandardTensors<dim>::I +
>> (alpham + (alphas - alpham) * phi) * (theta1 - thetaeq1) *
>> StandardTensors<dim>::I;
>>
>> // Transformation Strain Tensor
>> const SymmetricTensor<2, dim> epsTransformation = (1.0/3.0) * 0.02
>> * StandardTensors<dim>::I *
>> (1 - phi);
>> //The Initial Strain
>> const SymmetricTensor<2, dim> epsInitial = SymmetricTensor<2,
>> dim>(epsThermal) +
>> (1-phi) * 0.02 * StandardTensors<dim>::I;
>>
>> // Elastic Strain Tensor
>> const SymmetricTensor<2, dim> epsElastic = symmetrize(Tensor<2,
>> dim>(eps)-
>> SymmetricTensor<2,
>> dim>(epsTransformation)-
>> SymmetricTensor<2,
>> dim>(epsInitial) -
>> SymmetricTensor<2, dim>(epsThermal));
>> // The Deviatoric Strain
>> const SymmetricTensor<2, dim> epsDeviatoric = symmetrize(Tensor<2,
>> dim>(eps)-(1.0/3.0) * eps_V*
>> StandardTensors<dim>::I);
>>
>> material->update_material_data(F, eps, phi, dphi,
>> doubleWellBarrier, etaBVG, etaBVM
>> );
>> F_inv = invert(F);
>> eps_V = material->get_eps_V();
>> elasStress = material->get_elasStress();
>> // surfTenStress = material->get_surfTenStress(); // It is
>> exactly equal to initial stress...
>> totalStress = material->get_totalStress(); // Extracting
>> total Cauchy stress
>> Jc = material->get_Jc();
>> // Extracting elasticity stiffnes tensor
>> driving_force_noStress = material->get_driving_force_noStress();
>> // Extracting driving force with no stress
>>
>> eta_boundary_valueG = material->get_eta_boundary_valueG();
>> // Extracting the boundary value for G-L
>> equation
>> eta_boundary_valueM =
>> material->get_eta_boundary_valueM();
>> // Extracting the boundary value for mechanics
>>
>> const double temp = -(1.0/(1+eps_V))*0.02*mean_elas_stress;
>> const double temp1 = temp + 3*(1.0/(1+eps_V))*mean_elas_stress*
>> (alphas - alpham)*(theta1 - thetaeq1);
>> const double temp2 = temp1 -
>> ((1.0/(1+eps_V))*(0.5*(lambda1-lambda0)*
>> SymmetricTensor<2, dim>(epsElastic) *
>> SymmetricTensor<2, dim>(epsElastic)));
>> const double temp3 = -(1.0/(1+eps_V))*mu*SymmetricTensor<2,
>> dim>(epsDeviatoric)*
>> SymmetricTensor<2, dim>(epsDeviatoric);
>>
>> // Computing total driving force
>> get_driv_forc = temp2 * temp3 * dphi + driving_force_noStress;
>>
>> // computing the total boundary value for the G-L equation
>> get_eta_bou_valG = eta_boundary_valueG;
>>
>> // computing the total boundary value for the Mechanics
>> get_eta_bou_valM = eta_boundary_valueM;
>> Assert(determinant(F_inv) > 0, ExcInternalError());
>> }
>>
>> double get_det_F() const
>> {
>> return material->get_det_F();
>> }
>>
>> double get_eps_V() const
>> {
>> return material->get_eps_V();
>> }
>>
>> const Tensor<2, dim> &get_F_inv() const
>> {
>> return F_inv;
>> }
>>
>> const SymmetricTensor<2, dim> &get_elasStress() const
>> {
>> return elasStress;
>> }
>>
>> const SymmetricTensor<2, dim> &get_totalStress() const
>> {
>> return totalStress;
>> }
>>
>> // const SymmetricTensor<2, dim> &get_surfTenStress() const
>> // {
>> // return surfTenStress;
>> // }
>>
>> const SymmetricTensor<4, dim> &get_Jc() const
>> {
>> return Jc;
>> }
>>
>> double get_driving_force() const
>> {
>> return get_driv_forc;
>> }
>>
>> double get_meanElStress() const
>> {
>>
>> return mean_elas_stress;
>> }
>>
>> // The boundary value for the G-L equation
>> double get_eta_orderG() const
>> {
>> return get_eta_bou_valG;
>> }
>>
>> // The boundary value for the Mechanics
>> double get_eta_orderM() const
>> {
>> return get_eta_bou_valM;
>> }
>> private:
>> Material_Variables<dim> *material;
>> Tensor<2, dim> eps;
>> Tensor<2, dim> F_inv;
>> SymmetricTensor<2, dim> epsTransformation;
>> SymmetricTensor<2, dim> epsThermal;
>> SymmetricTensor<2, dim> epsInitial;
>> SymmetricTensor<2, dim> epsElastic;
>> SymmetricTensor<2, dim> epsDeviatoric;
>> double eps_V;
>> Tensor<2, dim> elasStress;
>>
>> Tensor<2, dim> totalStress;
>> // SymmetricTensor<2, dim> initialStress;
>> SymmetricTensor<4, dim> Jc;
>>
>
--
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) 2010 - 2015 by the deal.II authors and
* & Jean-Paul Pelteret and Andrew McBride
*
* 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.
*
* ---------------------------------------------------------------------
*
* Authors: Jean-Paul Pelteret, University of Cape Town,
* Andrew McBride, University of Erlangen-Nuremberg, 2010
*/
#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/dofs/dof_handler.h> // I add this....
#include <deal.II/dofs/dof_accessor.h> // I add this...
#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/grid/manifold_lib.h> // I add this...
#include <deal.II/grid/tria_accessor.h> // I add this...
#include <deal.II/grid/tria_iterator.h> // I add this...
#include <deal.II/fe/fe_dgp_monomial.h>
#include <deal.II/fe/fe_values_extractors.h> // I add this by my self...
#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/compressed_sparsity_pattern.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 <iostream>
#include <fstream>
#include <cmath>
//////////////////////////////////////////////////////////////////
namespace melting
{
using namespace dealii;
// INPUT OF PARAMETERS FOR AL NANOPARTICLES
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_edge;
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 edge", "32",
Patterns::Integer(0),
"Number of elements per long edge of the disc");
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_edge = prm.get_integer("Elements per edge");
scale = prm.get_double("Grid scale");
}
prm.leave_subsection();
}
/////////////////////////////////////////////////
struct Materials
{
double rout; // The external radius of disk
double lambda0; // Bulk modulus of liquid phase
double mu0; // Shear modulus of liquid phase
double lambda1; // Bulk modulus of solid phase
double mu1; // Shear modulus of solid phase
double chi; // Kinetic coefficient, Positive Definite
double betta; // Gradient energy coefficient, Positive Definite
double H; // Transformation heat
double thetain; // Initial temperature
double thetaeq; // Equilibrium temperature
double thetacr; // Critical temperature: 0.8*thetaeq
double theta; // Temperature of the particle
double pout; // Atmospheric pressure
double alphas; // Solid thermal expansion coefficient
double alpham; // Liquid thermal expansion coefficient
double etmr; // Transformation stress of melting in r-direction
double gammasv; // Surface energy between solid and vapor
double delg; // The difference between surface tension between solid and liquid
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("The external radius of disk", "30e-9", /* */
Patterns::Double(),
"The external radius of disk");
prm.declare_entry("Bulk modulus of liquid phase", "41.3e9",
Patterns::Double(0.0),
"Bulk modulus of liquid phase");
prm.declare_entry("Shear modulus of liquid phase", "0.0",
Patterns::Double(),
"Shear 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(0.0),
"Shear modulus of solid phase");
prm.declare_entry("Kinetic coefficient", "532.0",
Patterns::Double(0.0),
"Kinetic coefficient");
prm.declare_entry("Gradient energy coefficient", "3.21e-10",
Patterns::Double(0.0),
"Gradient energy coefficient");
prm.declare_entry("Transformation heat", "933.57e6",
Patterns::Double(0.0),
"Transformation heat");
prm.declare_entry("Initial temperature", "293.15",
Patterns::Double(0.0),
"Initial temperature");
prm.declare_entry("Equilibrium temperature", "933.67",
Patterns::Double(0.0),
"Equilibrium temperature");
prm.declare_entry("Critical temperature", "746.94",
Patterns::Double(0.0),
"Critical temperature");
prm.declare_entry("Temperature of the particle", "930.0",
Patterns::Double(0.0),
"Temperature of the particle");
prm.declare_entry("Atmospheric pressure", "0.0",
Patterns::Double(),
"Atmospheric pressure");
prm.declare_entry("Solid thermal expansion coefficient", "3.032e-5",
Patterns::Double(0.0),
"Solid thermal expansion coefficient");
prm.declare_entry("Liquid thermal expansion coefficient", "4.268e-5",
Patterns::Double(0.0),
"Liquid thermal expansion coefficient");
prm.declare_entry("Transformation stress of melting in r-direction", "0.02",
Patterns::Double(0.0),
"Transformation stress of melting in r-direction");
prm.declare_entry("Surface energy between solid and vapor", "1.05",
Patterns::Double(0.0),
"Surface energy between solid and vapor");
prm.declare_entry("The difference between surface tension between solid and liquid", "0.1292",
Patterns::Double(0.0),
"The difference between surface tension between solid and liquid");
}
prm.leave_subsection();
}
void Materials::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Material properties");
{
rout = prm.get_double("The external radius of disk");
lambda0 = prm.get_double("Bulk modulus of liquid phase");
mu0 = prm.get_double("Shear modulus of liquid phase");
lambda1 = prm.get_double("Bulk modulus of solid phase");
mu1 = prm.get_double("Shear modulus of solid phase");
chi = prm.get_double("Kinetic coefficient");
betta = prm.get_double("Gradient energy coefficient");
H = prm.get_double("Transformation heat");
thetain = prm.get_double("Initial temperature");
thetaeq = prm.get_double("Equilibrium temperature");
thetacr = prm.get_double("Critical temperature");
theta = prm.get_double("Temperature of the particle");
pout = prm.get_double("Atmospheric pressure");
alphas = prm.get_double("Solid thermal expansion coefficient");
alpham = prm.get_double("Liquid thermal expansion coefficient");
etmr = prm.get_double("Transformation stress of melting in r-direction");
gammasv = prm.get_double("Surface energy between solid and vapor");
delg = prm.get_double("The difference between surface tension between solid and liquid");
}
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 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", "2.5e-12",
Patterns::Double(),
"End time");
prm.declare_entry("Time step size", "0.5e-14",
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 Tensor
template <int dim>
class StandardTensors
{
public:
static const SymmetricTensor<2, dim> I;
};
// Static Constant SymmetricTensor<2, dim> for defining various strain;
template <int dim>
const SymmetricTensor<2, dim>
StandardTensors<dim>::I = unit_symmetric_tensor<dim>();
// Define time step, current time step, 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)
{}
virtual ~Time()
{}
double current() const
{
return time_current;
}
double end() const
{
return time_end;
}
double get_delta_t() const
{
return delta_t;
}
unsigned int get_timestep() const
{
return timestep;
}
void increment()
{
time_current += delta_t;
++timestep;
}
private:
unsigned int timestep;
double time_current;
const double time_end;
const double delta_t;
};
/////////////////////////////////////////////////////////////////////
template<int dim>
class Material_Variables
{
public:
Material_Variables(const double lambda0,
const double mu0,
const double lambda1,
const double mu1,
const double chi,
const double H,
// const double betta,
const double thetain,
const double thetaeq,
const double thetacr,
const double theta,
const double gammasv)
:
det_F(1.0),
eps_V(0.0), // For melting, it is equal to J=1+epsVol
eps(Tensor<2, dim>()), // 0.5*[grad_u + grad_u^T]
// epsVol(0.0), // epsVol = trace(eps)-(double contraction between total strain and unit tensor)
// epsAvg(0.0) // The average of strain; epsAvg = epsVol/3
// epsDeviatoric(SymmetricTensor<2, dim>()), // The deviatoric part of the strain tensor
// epsElastic(SymmetricTensor<2, dim>()), // The elastic strain
// epsInitial(SymmetricTensor<2, dim>()), // The initial strain
lambda00(lambda0),
mu00(mu0),
lambda11(lambda1),
mu11(mu1),
chi1(chi),
H1(H),
// betta1(betta),
thetain1(thetain),
thetaeq1(thetaeq),
thetacr1(thetacr),
theta1(theta),
gammasv1(gammasv)
{}
~Material_Variables()
{}
void update_material_data
( const Tensor<2, dim> &F, const Tensor<2, dim> &eps, const double phi, const double dphi,
const double doubleWellBarrier, const double etaBVG, const double etaBVM
)
{
det_F = determinant(F);
eps_V = trace(eps); // determinant of total F- The ratio of densities
// eps = 0.5*(Grad_u + transpose(Grad_u)); // eps = 1/2(grad_u + grad_u^T)
// epsVol = trace(eps); // Volumetric strain
// epsAvg = epsvol/3.0; // The mean of strain
// epsDev = eps - epsAvg; // The deviatoric part of the strain tensor
// epsElastic = epsVol - 3 * epsThermal; // Elastic Strain
// det_F = 1+epsVol; // The determinant of the deformation tensor(J)
mu = mu11*phi; // The Lame constant related to shear stress()
lambda = lambda00 + (lambda11 - lambda00)*phi; // The Lame constant related to both shear and Bulk modulus
A = 3*H1*(1-(thetacr1/thetaeq1)); // Compute barrier height
gam1 = (gammasv1 - delg1) + delg1 * phi; // The Surface Tension between Solid and Vapor
doubleWellBarrier1 = doubleWellBarrier; // The Helmholtz free energy according to double well barrier
etaBVG1 = etaBVG; // The boundary condition on the curved edge for G-L equation
etaBVM1 = etaBVM; // The boundary condition on the curved edge for mechanic problem
// thermalFreeEnergy1 = thermalFreeEnergy; // The thermal contribution to the Helmholtz free energy
dphi1 = dphi; // d_phi/d_eta
// Assert(eps_V > 0, ExcInternalError());
}
// The total Cauchy stress tensor in my problem is decomposed into
// elastic stress tensor and surface tension at interfaces.
// First we define the elastic stress tensor, then the surface tension stress
// at interfaces. The submission of these stress tensor form the total stress
// tensor.
SymmetricTensor<2, dim> get_elasStress() // The Elastic stress
{
return get_elastic_stress();
}
// SymmetricTensor<2, dim> get_meanElasticStress() // The mean elastic stress
// {
//
// return get_mean_elastic_stress();
// }
// SymmetricTensor<2, dim> get_surfTenStress() // The Surface_Tension stress at interfaces
// { // It is exactly equal to initial stress.
// return get_surface_tension_stress();
// }
SymmetricTensor<2, dim> get_totalStress() // The total Cauchy stress tensor
{
return get_total_stress();
}
// SymmetricTensor<2, dim> get_initialStress() // The initial stress for mechanics problem- Surface Tension Stress
// {
// return get_initial_stress();
// }
SymmetricTensor<4, dim> get_Jc() const // Comupte the Fourth Order Modulus Tensor
{
return get_Jc_modulus();
}
double get_det_F() const
{
return det_F;
}
double get_eps_V() const
{
return eps_V;
}
double get_meanElasticStress() const // It should be noted that the mean elastic stress is a scalar, not a tensor...
{
return get_mean_elastic_stress ();
}
double get_driving_force_noStress() const // Driving force without mechanical part (stress)
{
return get_driv_forc_noStress ();
}
double get_eta_boundary_valueG() const // Drive the boundary value for the order parameter
{
return get_bound_valueG ();
}
double get_eta_boundary_valueM () const // Drive the boundary value for mechanics
{
return get_bound_valueM ();
}
protected:
double det_F;
double eps_V;
Tensor<2, dim> eps;
Tensor<2, dim> epsThermal;
Tensor<2, dim> epsDeviatoric; // I add this by my self...
Tensor<2, dim> epsInitial; // I add this by my self...
Tensor<2, dim> epsElastic; // I add this by my self...
// double eta; // I add this parameter here.
// Tensor<1, dim> Grad_eta; // I add this function here.
double lambda00;
double mu00;
double lambda11;
double mu11;
double chi1;
double H1;
double betta1;
double delg1;
double thetain1;
double thetaeq1;
double thetacr1;
double theta1;
double gammasv1;
double lambda;
double mu;
double A;
double gam1;
double dphi1;
double doubleWellBarrier1;
double etaBVG1;
double etaBVM1;
double alphas;
double alpham;
double A1;
double phi;
// Compute the elastic stress tensor
SymmetricTensor<2, dim> get_elastic_stress() const
{
SymmetricTensor<2, dim> elastic_stress = symmetrize(0.04 * lambda*
StandardTensors<dim>::I +
2*mu*phi*Tensor<2, dim>(epsDeviatoric));
return elastic_stress;
}
// Compute the surface-tension stress at interfaces. First I define a 2rank symmetric tensor for the Dyadic product of grad_eta.
// SymmetricTensor<2, dim> get_surface_tension_stress() const // It is initial stress...
// {
// const Tensor<1, dim> Grad_eta;
// const SymmetricTensor<2, dim> grad_eta_x_grad_eta = outer_product(Grad_eta, Grad_eta);
// SymmetricTensor<2, dim> surface_tension_stress = 0.5*betta1*SymmetricTensor<2, dim>(grad_eta_x_grad_eta) +
// A*std::pow(eta,2)*std::pow(1-eta,2)*
// (Tensor<2, dim> (StandardTensors<dim>::I)-betta1 * SymmetricTensor<2, dim>(grad_eta_x_grad_eta));
// return surface_tension_stress;
// }
// Compute the total stress tensor
SymmetricTensor<2, dim> get_total_stress() const
{
const SymmetricTensor<2, dim> elastic_stress = get_elastic_stress();
// SymmetricTensor<2, dim> surface_tension_stress = get_surface_tension_stress();
SymmetricTensor<2, dim> total_stress = SymmetricTensor<2, dim>(elastic_stress);
return total_stress;
}
// Compute the elasticity stiffness 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;
}
double get_mean_elastic_stress() const
{
const SymmetricTensor<2, dim> elastic_stress = get_elastic_stress();
const double mean_elastic_stress = StandardTensors<dim>::I
* SymmetricTensor<2, dim>(elastic_stress) / 3.0;
return mean_elastic_stress;
}
// Computing the driving force without mechanical effects
double get_driv_forc_noStress() const
{
return -A*doubleWellBarrier1-H1*((theta1/thetaeq1) - 1)*dphi1;
}
double get_bound_valueG () const
{
return etaBVG1; // getting the boundary value for the G-L equation
}
double get_bound_valueM () const
{
return etaBVM1; // getting the boundary value for the mechanics
}
};
// Updates the Quadrature Point History
template <int dim>
class PointHistory
{
public:
PointHistory()
:
material(NULL),
F_inv(StandardTensors<dim>::I),
elasStress(SymmetricTensor<2, dim>()),
// surfTenStress(SymmetricTensor<2, dim>()),
totalStress(SymmetricTensor<2, dim>()),
Jc(SymmetricTensor<4, dim>())
{}
virtual ~PointHistory()
{
delete material;
material = NULL;
}
void setup_lqp (const Parameters::AllParameters ¶meters)
{
material = new Material_Variables<dim>(parameters.lambda0,
parameters.mu0, parameters.lambda1, parameters.mu1, parameters.chi, parameters.H,
parameters.thetain, parameters.thetacr, parameters.thetaeq, parameters.theta, parameters.gammasv);
update_values(Tensor<2, dim>(), double ());
//Tensor<1, dim>());
}
void update_values (const Tensor<2, dim> &Grad_u_n, const double eta
// const Tensor<1, dim> &Grad_eta
)
{
const Tensor<2, dim> F = (Tensor<2, dim>(StandardTensors<dim>::I) + Grad_u_n); // total F
const Tensor<2, dim> eps = 0.5 * (Grad_u_n + transpose(Grad_u_n)); // Total Strain
// I add this line for using grad_eta in the problem...
// const SymmetricTensor<2, dim> grad_eta_x_grad_eta = outer_product(Grad_eta, Grad_eta);
const double phi = std::pow(eta, 2) * (3 - 2 * eta); // Function of order parameter
dphi = 6 * eta * (1 - eta); // Derivative of the interpolation function
// The derivative of double-well energy
const double doubleWellBarrier = 2*eta + 4 * std::pow(eta, 3) - 6*std::pow(eta, 2);
const double etaBVG = -6 * chi1 * delg1 * (1+eps_V) * dphi; // Boundary value for G-L equation
const double etaBVM = 2 * gam1 / 30e-9; // Boundary value for Mechanics
// Thermal Strain Tensor
const SymmetricTensor<2, dim> epsThermal = alphas * (thetaeq1 - thetain1) * StandardTensors<dim>::I +
(alpham + (alphas - alpham) * phi) * (theta1 - thetaeq1) * StandardTensors<dim>::I;
// Transformation Strain Tensor
const SymmetricTensor<2, dim> epsTransformation = (1.0/3.0) * 0.02 * StandardTensors<dim>::I *
(1 - phi);
//The Initial Strain
const SymmetricTensor<2, dim> epsInitial = SymmetricTensor<2, dim>(epsThermal) +
(1-phi) * 0.02 * StandardTensors<dim>::I;
// Elastic Strain Tensor
const SymmetricTensor<2, dim> epsElastic = symmetrize(Tensor<2, dim>(eps)-
SymmetricTensor<2, dim>(epsTransformation)-
SymmetricTensor<2, dim>(epsInitial) -
SymmetricTensor<2, dim>(epsThermal));
// The Deviatoric Strain
const SymmetricTensor<2, dim> epsDeviatoric = symmetrize(Tensor<2, dim>(eps)-(1.0/3.0) * eps_V*
StandardTensors<dim>::I);
material->update_material_data(F, eps, phi, dphi, doubleWellBarrier, etaBVG, etaBVM
);
F_inv = invert(F);
eps_V = material->get_eps_V();
elasStress = material->get_elasStress();
// surfTenStress = material->get_surfTenStress(); // It is exactly equal to initial stress...
totalStress = material->get_totalStress(); // Extracting total Cauchy stress
Jc = material->get_Jc(); // Extracting elasticity stiffnes tensor
driving_force_noStress = material->get_driving_force_noStress(); // Extracting driving force with no stress
eta_boundary_valueG = material->get_eta_boundary_valueG(); // Extracting the boundary value for G-L equation
eta_boundary_valueM = material->get_eta_boundary_valueM(); // Extracting the boundary value for mechanics
const double temp = -(1.0/(1+eps_V))*0.02*mean_elas_stress;
const double temp1 = temp + 3*(1.0/(1+eps_V))*mean_elas_stress*
(alphas - alpham)*(theta1 - thetaeq1);
const double temp2 = temp1 - ((1.0/(1+eps_V))*(0.5*(lambda1-lambda0)*
SymmetricTensor<2, dim>(epsElastic) * SymmetricTensor<2, dim>(epsElastic)));
const double temp3 = -(1.0/(1+eps_V))*mu*SymmetricTensor<2, dim>(epsDeviatoric)*
SymmetricTensor<2, dim>(epsDeviatoric);
// Computing total driving force
get_driv_forc = temp2 * temp3 * dphi + driving_force_noStress;
// computing the total boundary value for the G-L equation
get_eta_bou_valG = eta_boundary_valueG;
// computing the total boundary value for the Mechanics
get_eta_bou_valM = eta_boundary_valueM;
Assert(determinant(F_inv) > 0, ExcInternalError());
}
double get_det_F() const
{
return material->get_det_F();
}
double get_eps_V() const
{
return material->get_eps_V();
}
const Tensor<2, dim> &get_F_inv() const
{
return F_inv;
}
const SymmetricTensor<2, dim> &get_elasStress() const
{
return elasStress;
}
const SymmetricTensor<2, dim> &get_totalStress() const
{
return totalStress;
}
// const SymmetricTensor<2, dim> &get_surfTenStress() const
// {
// return surfTenStress;
// }
const SymmetricTensor<4, dim> &get_Jc() const
{
return Jc;
}
double get_driving_force() const
{
return get_driv_forc;
}
double get_meanElStress() const
{
return mean_elas_stress;
}
// The boundary value for the G-L equation
double get_eta_orderG() const
{
return get_eta_bou_valG;
}
// The boundary value for the Mechanics
double get_eta_orderM() const
{
return get_eta_bou_valM;
}
private:
Material_Variables<dim> *material;
Tensor<2, dim> eps;
Tensor<2, dim> F_inv;
Tensor<2, dim> epsTransformation;
Tensor<2, dim> epsThermal;
Tensor<2, dim> epsInitial;
Tensor<2, dim> epsElastic;
Tensor<2, dim> epsDeviatoric;
double eps_V;
Tensor<2, dim> elasStress;
// SymmetricTensor<2, dim> surfTenStress; // It is exactly the same with initial stress...
Tensor<2, dim> totalStress;
// SymmetricTensor<2, dim> initialStress;
SymmetricTensor<4, dim> Jc;
// Tensor<2, dim> Fe;
double thetain1;
double thetacr1;
double thetaeq1;
double theta1;
double lambda0;
double lambda1;
double chi1;
double H1;
double betta1;
double alphas;
double alpham;
double mu;
double delg1;
double gam1;
double driving_force_noStress;
double mean_elas_stress;
double eta_boundary_valueG; // the boundary value for the G-L equation
double eta_boundary_valueM; // the boundary value for the Mechanics
double dphi;
double get_driv_forc;
double get_eta_bou_valG; // Call the boundary value for the G-L equation
double get_eta_bou_valM; // Call the boundary value for the Mechanics
};
// The Mechanical part of the problem
template <int spacedim>
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_mass;
struct ScratchData_mass;
struct PerTaskData_laplace;
struct ScratchData_laplace;
struct PerTaskData_PF_RHS;
struct ScratchData_PF_RHS; */
struct PerTaskData_UQPH;
struct ScratchData_UQPH;
static const unsigned int dim = spacedim - 1; // I add this for computing on curved domain...
// Function creates geometry and mesh
void
make_grid();
// Sets of all the vectors and matrices for Ginzburg-Landau eqn and the mechanics problem
void
system_setup();
// Assemble the tangent stiffness for N-R problem for solving nodal displacements
void
assemble_system_tangent();
void
assemble_system_tangent_one_cell(const typename
DoFHandler<dim, spacedim>::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 or the residual
void
assemble_system_rhs();
void
assemble_system_rhs_one_cell(const typename
DoFHandler<dim, spacedim>::active_cell_iterator &cell,
ScratchData_RHS &scratch,
PerTaskData_RHS &data);
void
copy_local_to_global_rhs(const PerTaskData_RHS &data);
/*
void
assemble_system_mass();
void
assemble_system_mass_one_cell(const typename
DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_mass &scratch,
PerTaskData_mass &data);
void
copy_local_to_global_mass(const PerTaskData_mass &data);
void
assemble_system_laplace();
void
assemble_system_laplace_one_cell(const typename
DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_laplace &scratch,
PerTaskData_laplace &data);
void
copy_local_to_global_laplace(const PerTaskData_laplace &data);
void
assemble_system_PF_rhs();
void
assemble_system_PF_rhs_one_cell(const typename
DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_PF_RHS &scratch,
PerTaskData_PF_RHS &data);
void
copy_local_to_global_PF_rhs(const PerTaskData_PF_RHS &data);
*/
// Imposes displacement constraint
void
make_constraints(const int &it_nr);
// Sets up quadrature point history
void
setup_qph();
// Functions for updating the quadrature point history
void update_qph_incremental(const BlockVector<double> &solution_delta);
//const Vector<double> &solution_delta_eta); // I add one member function
void
update_qph_incremental_one_cell(const typename
DoFHandler<dim, spacedim>::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 updated values of the order parameter
void
solve_eta(Vector<double> &solution_delta_eta);
// Computes the linear problem for obtaining total order parameter(I add this member function by my self....)
// void
// get_total_solution_eta(Vector<double> &solution_delta_eta_total);
// Computes updated values of the gradient of the order parameter
// void
// solve_grad_eta(Tensor<2, dim> &solution_grad_delta_eta); // I should compute the gradient of the order parameter....
// Computes the right hand side in the order parameter equation
void
compute_eta_rhs (Vector<double> &force_eta) const;
void
output_results() const;
Parameters::AllParameters parameters;
double vol_reference; // Volume of the reference config
double vol_current; // Volume of the current config
// I add two more triangulation for defining my geometry in two dimension...
Triangulation<dim, spacedim> 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, spacedim> fe; // Fe object
DoFHandler<dim, spacedim> dof_handler_ref; // We have used two dof_handler: one for mechanics prob
// and the other for Ginzburg-Landau order parameter
MappingQ<dim, spacedim> mapping;
const unsigned int dofs_per_cell; // No of dofs per cell for the mechanics problem
// unsigned int n_u = 0; // I add this line by my self for displacement along horizontal axis....
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 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(order parameter)
FE_Q<dim, spacedim> fe_eta; // Fe object for eta
DoFHandler<dim, spacedim> dof_handler_eta; // Another dof_handler for eta
const unsigned int dofs_per_cell_eta; // Dof per eta cell
SparsityPattern sparsity_pattern_eta;
SparseMatrix<double> mass_matrix; // Mass matrix out of Ginzburg-Landau eqn
SparseMatrix<double> laplace_matrix; // Laplace matrix out of Ginzburg-Landau eqn
Vector<double> solution_eta; // Solution vector for eta
// Tensor<2, dim> solution_grad_eta; // Solution vector for grad_eta+++++++++++++++++
Vector<double> tmp_vector; // A vector used for solving 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 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
{
/*
if ((p[0] > 2.5e-9) && (p[1] > 2.5e-9) && (p[0] < 7.5e-9) && (p[1] < 7.5e-9))
return 1;
else
return 0;
return
-std::max(std::min(-tanh((sqrt((p[0]-3.0e-9)*(p[0]-3.0e-9)+(p[1]-3.0e-9)*(p[1]-3.0e-9))-1.0e-9)/0.4e-9),0.0),
std::min(-tanh((sqrt((p[0]-7.0e-9)*(p[0]-7.0e-9)+(p[1]-3.0e-9)*(p[1]-3.0e-9))-1.0e-9)/0.4e-9),0.0));
*/
return 0;
// return std::max(tanh((3.5e-9-sqrt((p(0)-4.0e-9)*(p(0)-4.0e-9)+(p(1)-4.0e-9)*
// (p(1)-4.0e-9)))/0.4e-9),0.0);
//return (0.5*std::tanh(3.0*(p[0]-5.0e-9)/7.07e-10)+0.5);
}
///////////////////////////////////////////////////
// Constructor
template <int spacedim>
Solid<spacedim>::Solid(const std::string &input_file)
:
parameters(input_file),
triangulation(Triangulation<dim, spacedim>::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, spacedim>(parameters.poly_degree), dim), // Displacement
dof_handler_ref(triangulation),
mapping(parameters.poly_degree), // Using mapping for curved domain
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), // Order Parameter
dof_handler_eta (triangulation),
dofs_per_cell_eta(fe_eta.dofs_per_cell)
{
// determine_component_extractor();
}
// Destructor
template <int spacedim>
Solid<spacedim>::~Solid()
{
dof_handler_ref.clear(); // Dof Handler for displacement
dof_handler_eta.clear(); // Dof Handler for order parameter
}
//////////////////////////////////////////////////
template <int spacedim>
void Solid<spacedim>::run()
{
make_grid(); // 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> force_eta (solution_eta.size());
// 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;
// solution_grad_delta_eta = 0.0; // I add this by my self...
force_eta = 0.0;
tmp_vector = 0.0;
// Vector<double> tmp_vector(solution_eta.size());
/*
MatrixCreator::create_mass_matrix(dof_handler_eta,
QGauss<dim>(fe_eta.degree+1), //change here deg
mass_matrix);
MatrixCreator::create_laplace_matrix(dof_handler_eta,
QGauss<dim>(fe_eta.degree+1), //change here deg
laplace_matrix);
*/
laplace_matrix.vmult (tmp_vector, solution_eta);
tmp_vector *= -parameters.betta;
compute_eta_rhs (force_eta);
std::cout << force_eta.l2_norm () << std::endl;
std::cout << tmp_vector.l2_norm () << std::endl;
tmp_vector += force_eta;
tmp_vector *= parameters.chi;
tmp_vector *= parameters.delta_t;
solve_eta (solution_delta_eta); // Solving the Order parameter
/* Vector<double> junk (solution_eta.size());
mass_matrix.vmult (junk, solution_delta_eta);
std::cout << junk.l2_norm () << std::endl;
*/
solution_eta += solution_delta_eta; // Total eta
update_qph_incremental(solution_delta);
//solution_delta_eta); // Updating the quadrature point history. Add one memeber
// 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 spacedim>
struct Solid<spacedim>::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 spacedim>
struct Solid<spacedim>::ScratchData_K
{
FEValues<dim, spacedim> 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, spacedim> &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 spacedim>
struct Solid<spacedim>::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 spacedim>
struct Solid<spacedim>::ScratchData_RHS
{
FEValues<dim, spacedim> fe_values_ref;
FEFaceValues<dim, spacedim> fe_face_values_ref;
std::vector<std::vector<double> > Nx;
std::vector<std::vector<SymmetricTensor<2, dim> > > symm_grad_Nx;
// Here I define some object for implying boundary condition of G-L equation...
ScratchData_RHS(const FiniteElement<dim, spacedim> &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;
}
}
}
};
// Next two functions are used in updating the quadrature point history with new stress, driving force
// and other quantities such as eta boundary values for mechanics and G-L equation...
/////////////////////////////////////////////////////
template <int spacedim>
struct Solid<spacedim>::PerTaskData_UQPH
{
void reset()
{}
};
////////////////////////////////////////////////////
template <int spacedim>
struct Solid<spacedim>::ScratchData_UQPH
{
const BlockVector<double> &solution_total;
// const std::vector<double> &solution_eta;
std::vector<Tensor<2, dim> > solution_grads_u_total;
std::vector<double> solution_eta_total; // Check if required
// std::vector<Tensor<2, dim> > solution_grads_eta_total; // I add this by my self......
FEValues<dim, spacedim> fe_values_ref;
FEValues<dim, spacedim> fe_values_eta;
ScratchData_UQPH(const FiniteElement<dim, spacedim> &fe_cell,
const UpdateFlags uf_cell,
const FiniteElement<dim, spacedim> &fe_cell_eta,
const UpdateFlags uf_cell_eta,
const QGauss<dim> &qf_cell,
const BlockVector<double> &solution_total)
// std::vector<double> solution_eta_total) // I add this by my self...
:
solution_total(solution_total),
solution_grads_u_total(qf_cell.size()),
// solution_eta(solution_eta),
solution_eta_total(qf_cell.size()),
// solution_grads_eta_total(qf_cell.size()), // I add this line by my self....
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(rhs.solution_eta),
solution_eta_total(rhs.solution_eta_total),
// solution_grads_eta_total(rhs.solution_grads_eta_total), // I add this line by my self...
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;
// solution_grads_eta_total[q] = 0.0; // I add this line by my self.......
}
}
};
/////////////////////////////////////////////////////////////
// This function computes the right hand terms (non-linear) in the Ginzburg-Landau equation
template <int spacedim>
void Solid<spacedim>::compute_eta_rhs (Vector<double> &nl_term) const
{
nl_term = 0; // Initialize the global non-linear term in Eq. 1
// FEFaceValues<dim, spacedim> fe_face_values_eta; // I add this for computing order parameter for each face
const QGauss<dim> quadrature_formula_eta (fe_eta.degree+1); // Change here
FEValues<dim, spacedim> fe_values_eta (mapping, 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();
Vector<double> local_nl_term (dofs_per_cell_eta); // Define the non-linear term in Eq. 1 for each element
std::vector<types::global_dof_index> local_dof_indices_eta (dofs_per_cell_eta);
// std::vector<double> old_data_values (n_q_points_eta);
typename DoFHandler<dim, spacedim>::active_cell_iterator
cell_eta = dof_handler_eta.begin_active(),
endc = dof_handler_eta.end();
// Loop to compute the non-linear term in the right hand side for each element
for (; cell_eta!=endc; ++cell_eta)
{
local_nl_term = 0; // Initialize the non-linear term in Eq. 1 for an element
fe_values_eta.reinit (cell_eta);
PointHistory<dim> *lqph =
reinterpret_cast<PointHistory<dim>*>(cell_eta->user_pointer());
for (unsigned int q_point=0; q_point<n_q_points_eta; ++q_point)
{ // Add the three following variable for the cylindrical coordinate system
double pii = numbers::PI;
double r_value = fe_values_eta.quadrature_point(q_point)[0];
const double driving_force = lqph[q_point].get_driving_force();
for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
{ // The term for order parameter boundary value is considered in this function in
// Cylindrical coordinate system.......
local_nl_term(i) += 2 * pii * r_value * driving_force *
fe_values_eta.shape_value(i, q_point) *
fe_values_eta.JxW(q_point);
}
}
cell_eta->get_dof_indices (local_dof_indices_eta);
// Adding the nodal terms into the global non-linear vector
for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
nl_term(local_dof_indices_eta[i]) += local_nl_term(i);
}
// // Here, We imply the etaBoundaryValue for G-L equation, We first check to see it the cell face exists on a boundary
// // the etaBoundaryValue is applied and add the contribution if this is the case...
// for (unsigned int face_eta = 0; face_eta<GeometryInfo<dim>::faces_per_cell; ++face_eta)
// if (cell_eta->face(face_eta)->at_boundary() == true
// &&cell_eta->face(face_eta)->boundary_id() == 0)
// {
// const QGauss<dim> quadrature_formula_ref (fe.degree+1); // Change here
// PointHistory<dim> *lqph =
// reinterpret_cast<PointHistory<dim>*>(cell_eta->user_pointer());
// for (unsigned int q_point = 0; q_point<n_q_points_eta; ++q_point)
// { // We specify the etaBoundaryValue in the current configuration. For this problem, the etaBoundaryValue
// // is implied in cylindrical coordinate system on b_id0.
// const double eta_orderG = lqph[q_point].get_eta_orderG();
// const double r_value = fe_values_eta.quadrature_point(q_point)[0];
// const double eta_V = lqph[q_point].get_eta_V();
// for (unsigned int f_q_point = 0; f_q_point<n_q_points_f; ++f_q_point)
// { // Defining the pressure on the curved edge for mechanics...
// const double time_ramp = (time.current() / time.end());
// double pii = numbers::PI;
// const double pressure_on_boundaryG = eta_orderG * pii ;
//
// for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
// {
// const double Ei = fe_values_eta.shape_value(i, f_q_point);
// const double JxW = fe_values_eta.JxW(f_q_point);
// nl_term(i) -= 2 * pressure_on_boundaryG * r_value * Ei * JxW;
// }
//
// }
//
//
// }
//
// }
}
/////////////////////////////////////////////////////////////////////////////
template <int spacedim>
void Solid<spacedim>::make_grid()
{ // This is the new domain that I created for my problem...
static SphericalManifold<dim,spacedim> surface_description;
{
Triangulation<spacedim> volume_mesh;
GridGenerator::half_hyper_ball(volume_mesh);
std::set<types::boundary_id> boundary_ids;
boundary_ids.insert (0);
GridGenerator::extract_boundary_mesh (volume_mesh, triangulation,
boundary_ids);
}
triangulation.set_all_manifold_ids(0);
triangulation.set_manifold (0, surface_description);
triangulation.refine_global(4);
std::cout << "Surface mesh has " << triangulation.n_active_cells()
<< " cells."
<< std::endl;
vol_reference = GridTools::volume(triangulation);
vol_current = vol_reference;
std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
// First, I define one triangulation for defining one quarter_hyper_shell with the following dimension...
// Triangulation<dim> tria1;
// const Point<dim> center;
// const double inner_radius = 3e-9;
// const double outer_radius = 30e-9;
// GridGenerator::quarter_hyper_shell (tria1,
// center,
// inner_radius,
// outer_radius,
// 0,
// false);
// tria1.set_manifold (0);
//
// // Second, I define another triangulation for the second domain...
// Triangulation<dim> tria2;
// const Point<dim> center1;
// const double inner_radius1 = 0.0001e-9;
// const double outer_radius1 = 3e-9;
// GridGenerator::quarter_hyper_shell(tria2,
// center1,
// inner_radius1,
// outer_radius1,
// 0,
// false);
//
//// Finally, we merge two triangulation together for creating the whole domain.....
//
// GridGenerator::merge_triangulations (tria1, tria2, triangulation);
// triangulation.set_all_manifold_ids(0);
// const SphericalManifold<dim> manifold_description1(center1);
// triangulation.set_manifold (0, manifold_description1);
// GridTools::copy_boundary_to_manifold_id(triangulation);
// triangulation.refine_global (2);
// triangulation.set_manifold (0);
// for (typename Triangulation<dim> ::active_cell_iterator
// cell=triangulation.begin_active();
// cell!=triangulation.end(); ++cell)
// for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
// if (cell->face(f)->at_boundary())
// {
// const Point<dim> face_center = cell->face(f)->center();
// if (face_center[0] == 0)
// cell->face(f)->set_boundary_id (0); // faces on y-axis...
// else if (face_center[1] == 0)
// cell->face(f)->set_boundary_id (1); // faces on x-axis...
// else if (std::sqrt(face_center[0]*face_center[0]+
// face_center[1]*face_center[1]) == 1)
// cell->face(f)->set_boundary_id (2); // faces on the curved edge of the domain...
// }
// GridTools::scale(parameters.scale, triangulation);
}
///////////////////////////////////////////////////////////////////////////////////
template <int spacedim>
void Solid<spacedim>::system_setup()
{
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);
// I add the boundary condition for the x-axis here as no_normal_flux_constraints....
// std::set<types::boundary_id> 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);
DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block,
block_component);
// 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); // It is equal to n_u for the right hand side...
system_rhs.collect_sizes();
solution_n.reinit(dofs_per_block); // the initialization of n_u for displacement...
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();
{
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);
mass_matrix.reinit (sparsity_pattern_eta);
laplace_matrix.reinit (sparsity_pattern_eta);
MatrixCreator::create_mass_matrix(dof_handler_eta,
QGauss<dim>(fe_eta.degree+1), // Creates the mass matrix
mass_matrix);
MatrixCreator::create_laplace_matrix(dof_handler_eta,
QGauss<dim>(fe_eta.degree+1), // Creates the Laplace matrix
laplace_matrix);
solution_eta.reinit(dof_handler_eta.n_dofs());
// solution_grad_eta.reinit(dof_handler_eta.n_dofs()); // I add this by my self.....
tmp_vector.reinit(dof_handler_eta.n_dofs());;
}
// Sets up the quadrature point history
setup_qph();
timer.leave_subsection();
}
//////////////////////////////////////////////////////
template <int spacedim>
void Solid<spacedim>::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, spacedim>::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, spacedim>::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 spacedim>
void Solid<spacedim>::update_qph_incremental(const BlockVector<double> &solution_delta)
//const Vector<double> &solution_delta_eta) // I add this by my self...
// Check if solution_eta to be passed here
{
timer.enter_subsection("Update QPH data");
std::cout << " UQPH " << std::flush;
const BlockVector<double> solution_total(get_total_solution(solution_delta));
// std::vector<double> solution_eta_total(get_total_solution_eta(solution_delta_eta_total)); // I add this by my self...
const UpdateFlags uf_UQPH(update_gradients);
const UpdateFlags uf_UQPH_eta(update_values);
// update_gradients); // I should add this whenever i want to compute the gradient of eta...
PerTaskData_UQPH per_task_data_UQPH;
ScratchData_UQPH scratch_data_UQPH(fe, uf_UQPH, fe_eta, uf_UQPH_eta,
qf_cell, solution_total);
// solution_eta_total); // I add this by my self...
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 spacedim>
void
Solid<spacedim>::update_qph_incremental_one_cell(const typename
DoFHandler<dim, spacedim>::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_eta_total.size() == n_q_points,
ExcInternalError());
scratch.reset();
scratch.fe_values_ref.reinit(cell);
typename DoFHandler<dim, spacedim>::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, // Check if scratch is needed
scratch.solution_eta_total);
// scratch.fe_values_eta.get_function_gradients(scratch.solution_eta_total,
// scratch.solution_grads_eta_total); // I add this by my self....
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_eta_total[q_point]);
// scratch.solution_grads_eta_total[q_point]); // I add this line by my self...
}
// Solving the order parameter
template <int spacedim>
void
Solid<spacedim>::solve_eta (Vector<double> &solution_delta_eta)
{
SolverControl solver_control (1000, 1e-12*tmp_vector.l2_norm());
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(mass_matrix, 1.2);
cg.solve (mass_matrix, solution_delta_eta, tmp_vector,
preconditioner); // Should I compute the first derivative of order parameter here??
// return solver_control.last_step();
}
// Solves the N-R problem for displacement
template <int spacedim>
void
Solid<spacedim>::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;
// Const BlockVector<double> solution_total(get_total_solution(solution_delta)); // Added
update_qph_incremental(solution_delta);
// solution_delta_eta); // I add this by my self...
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 spacedim>
void Solid<spacedim>::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 spacedim>
void Solid<spacedim>::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 spacedim>
void Solid<spacedim>::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 spacedim>
void Solid<spacedim>::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 spacedim>
BlockVector<double>
Solid<spacedim>::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 spacedim>
void Solid<spacedim>::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 spacedim>
void Solid<spacedim>::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 spacedim>
void
Solid<spacedim>::assemble_system_tangent_one_cell(const typename
DoFHandler<dim, spacedim>::active_cell_iterator &cell,
ScratchData_K &scratch,
PerTaskData_K &data)
{
// const QGauss<dim> quadrature_formula_ref (fe.degree+1);
// FEValues<dim> fe_values_ref (fe, quadrature_formula_ref,
// update_values |
// update_JxW_values |
// update_quadrature_points |
// update_gradients); // I add this for computing displacement in cylindrical coordinate system
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();
// const double eps_V = lqph[q_point].get_eps_V(); // I add this by my self...
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)
{ // Add these variables for employing the cylindrical coordinate system to the equation....
double pii = numbers::PI;
// double z_value = fe_values_ref.quadrature_point(q_point)[1];
const Tensor<2, dim> totalStress = lqph[q_point].get_totalStress();
const SymmetricTensor<4, dim> Jc = lqph[q_point].get_Jc();
// const double det_F = lqph[q_point].get_det_F(); // I add this from anup code...
// const std::vector<double>
// &N = scratch.Nx[q_point];
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) += 2 * pii * symm_grad_Nx[i] * Jc
* symm_grad_Nx[j] * JxW;// The material contribution in cylindrical coordinate system
if (component_i == component_j) // Geometrical stress contribution in cylindrical coordinate system
data.cell_matrix(i, j) += 2 * pii * grad_Nx[i][component_i] * totalStress
* 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 spacedim>
void Solid<spacedim>::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);
// I add update_normal_vectors for computing the normal vector in the face for the boundary value related to the Mechanics
const UpdateFlags uf_face(update_values |
update_JxW_values
);
PerTaskData_RHS per_task_data(dofs_per_cell);
// I add two other scratch for boundary value of order parameter...
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 spacedim>
void Solid<spacedim>::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 spacedim>
void
Solid<spacedim>::assemble_system_rhs_one_cell(const typename
DoFHandler<dim, spacedim>::active_cell_iterator &cell,
ScratchData_RHS &scratch,
PerTaskData_RHS &data)
{
const QGauss<dim> quadrature_formula_ref (fe.degree+1);
FEValues<dim, spacedim> fe_values_ref (mapping, fe, quadrature_formula_ref,
update_values |
update_JxW_values |
update_quadrature_points |
update_gradients); // I add this for computing displacement in cylindrical coordinate system
// I add this for computing the boundary value of order parameter.....
// const QGauss<dim> quadrature_formula_eta (fe_eta.degree+1);
// FEValues<dim> fe_values_eta (fe_eta, quadrature_formula_eta,
// update_values |
// update_JxW_values |
// update_quadrature_points |
// update_gradients);
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)
// I make a change here...
scratch.symm_grad_Nx[q_point][k]
= 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 double det_F = lqph[q_point].get_det_F(); // I add this by my self...
double pii = numbers::PI;
double r_value = fe_values_ref.quadrature_point(q_point)[0];
// double z_value = fe_values_ref.quadrature_point(q_point)[1];
const SymmetricTensor<2, dim> totalStress = lqph[q_point].get_totalStress();
// const std::vector<double>
// &N = scratch.Nx[q_point];
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) -= 2 * pii * r_value * (symm_grad_Nx[i] * totalStress) * JxW;
else
Assert(i_group <= u_dof, ExcInternalError());
}
}
// Here, We imply the etaBoundaryValue for Mechanics, We first check to see it the cell face exists on a boundary
// the etaBoundaryValue is applied and add the contribution if this is the case...
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);
const QGauss<dim> quadrature_formula_ref (fe.degree+1); // Change here
for (unsigned int q_point = 0; q_point<n_q_points; ++q_point)
{ // We specify the etaBoundaryValue in the current configuration. For this problem, the etaBoundaryValue
// is implied in cylindrical coordinate system on b_id2.
// double pii = numbers::PI;
const double r_value = fe_values_ref.quadrature_point(q_point)[0];
const double eta_orderM = lqph[q_point].get_eta_orderM();
// const double eta_V = lqph[q_point].get_eta_V();
for (unsigned int f_q_point = 0; f_q_point<n_q_points_f; ++f_q_point)
{ // Defining the pressure on the curved edge for mechanics...
const double time_ramp = (time.current() / time.end());
double pii = numbers::PI;
const double pressure_on_boundaryM = eta_orderM * time_ramp * pii ;
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double Ei = scratch.fe_face_values_ref.shape_value(i, f_q_point);
const double JxW = scratch.fe_face_values_ref.JxW(f_q_point);
data.cell_rhs(i) -= 2 * pressure_on_boundaryM * r_value * Ei * JxW;
}
}
}
}
}
////////////////////////////////////////////////////////////////////////////
// Imposing the constraints
template <int spacedim>
void Solid<spacedim>::make_constraints(const int &it_nr)
{
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 z_displacement(1);
const FEValuesExtractors::Scalar phi_displacement(2);
{
const int boundary_id = 1;
if (apply_axialsymmetry_bc == true)
VectorTools::interpolate_boundary_values(mapping, dof_handler_ref,
boundary_id,
ZeroFunction<spacedim>(n_components),
constraints,
fe.component_mask(r_displacement)
);
// else
// VectorTools::interpolate_boundary_values(mapping, dof_handler_ref,
// boundary_id,
// ZeroFunction<dim>(n_components),
// constraints,
// fe.component_mask(r_displacement)
// );
}
// {
// const int boundary_id = 1;
// if (apply_axialsymmetry_bc == true)
// VectorTools::interpolate_boundary_values(dof_handler_ref,
// boundary_id,
// ZeroFunction<dim>(n_components),
// constraints,
// fe.component_mask(z_displacement)
// );
// else
// VectorTools::interpolate_boundary_values(dof_handler_ref,
// boundary_id,
// ZeroFunction<dim>(n_components),
// constraints,
// fe.component_mask(z_displacement)
// );
// }
// In this part, I fixed the center of the circle in r-direction and z-direction...
// {
// if (apply_axialsymmetry_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.0) &&
// (cell->vertex(vertex_index)(2)==0.5e-9))
// {
// std::cout << "r constr" << std::endl;
// constraints.add_line(cell->vertex_dof_index(vertex_index, 0));
// constraints.add_line(cell->vertex_dof_index(vertex_index, 2));
// }
// }
// }
// }
// }
constraints.close();
}
////////////////////////////////////////////////////////////////
template <int spacedim>
std::pair<unsigned int, double>
Solid<spacedim>::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 spacedim>
void Solid<spacedim>::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, DoFHandler<dim, spacedim> > data_out;
data_out.attach_dof_handler(dof_handler_ref);
data_out.add_data_vector(solution_n,
solution_name,
DataOut<dim,DoFHandler<dim, spacedim> >::type_dof_data,
data_component_interpretation);
// data_out.attach_dof_handler(dof_handler_eta); // I add this by my self...
data_out.add_data_vector (dof_handler_eta, solution_eta, "eta");
data_out.build_patches (mapping, mapping.get_degree()); // Check what is this- I add this by my self...
Vector<double> soln(solution_n.size());
for (unsigned int i = 0; i < soln.size(); ++i)
soln(i) = solution_n(i);
std::ostringstream filename;
filename << "solution-" << time.get_timestep() << ".vtk";
std::ofstream output(filename.str().c_str());
data_out.write_vtk(output);
}
}
/////////////////////////////////////////////////////////
// The main part of the program
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);
Solid<3> 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;
}
