Daniel, Let me explain the problem more. I heve parallelized a code based on step-40. The non-parallel code works for all the initial values and boundary conditions But the parallel one fails at the first call to solverCG with this error: "[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range" and the same story when I run the code on cluster.
For a couple of specific boundary conditions and initial values the problem resolves when I just replace the PreconditioerAMG with PreconditionerJacobi !! I feel the problem is with Petsc solver, that's why I want to examine if the code works with Trilinos. But I failed to implement trilinos as I mentioned previously. I am submitting my code along with a Parameter.prm file which includes all the parameters of the code. I appreciate if you could run the code on your machine and see if you can find out what's wrong with it. Thanks For example the following is the number of SolverCG iteration for the 20 Newoton-Raphson iterations after which On Sunday, October 16, 2016 at 10:39:49 AM UTC-5, Daniel Arndt wrote: > > In case there is an error with Trilinos, what is it? > > Am Sonntag, 16. Oktober 2016 17:27:26 UTC+2 schrieb Daniel Arndt: >> >> Hamed, >> >> >>>>> It still works with Petsc and Trilinos is not implemented at all. I >>>>> appreciate it if you could let me know what changes should be made in >>>>> step-40 to make it work with Trilinos. >>>>> >>>> What exactly, do you mean by that? So you have both PETSc and Trilinos >> installed and step-40 works. Are you saying that PETSc objects are used >> even if you use >> "using namespace ::LinearAlgebraTrilinos;" ? How did you find out? >> >> Best, >> Daniel >> > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
#include <deal.II/base/function.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/std_cxx11/shared_ptr.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/generic_linear_algebra.h>
namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
using namespace dealii::LinearAlgebraPETSc;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace dealii::LinearAlgebraTrilinos;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
}
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/compressed_sparsity_pattern.h>
#include <deal.II/lac/precondition_selector.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.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/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <iostream>
#include <fstream>
//////////////////////////////////////////////////////////////////
namespace PhaseField
{
using namespace dealii;
// INPUT OF PARAMETERS
namespace Parameters
{
struct FESystem
{
unsigned int poly_degree;
unsigned int quad_order;
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void FESystem::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Finite element system");
{
prm.declare_entry("Polynomial degree", "1",
Patterns::Integer(0),
"Displacement system polynomial order");
prm.declare_entry("Quadrature order", "2",
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 refinement;
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("Global refinement", "4",
Patterns::Integer(0),
"Global refinement level");
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");
{
refinement = prm.get_integer("Global refinement");
scale = prm.get_double("Grid scale");
}
prm.leave_subsection();
}
/////////////////////////////////////////////////
struct Materials
{
double lambda0; // austenite phase
double mu0; // austenite phase
double lambda1; // martensite phase
double mu1; // martensite phase
double L; // interface mobility
double beta; // gradient energy coefficient
double A0; // parameter for barrier height
double theta; // temperature
double thetac; // critical temperature
double thetae; // equilibrium temperature
double thetan; // Crank-Nicolson parameter
static void
declare_parameters(ParameterHandler &prm);
void
parse_parameters(ParameterHandler &prm);
};
void Materials::declare_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Material properties");
{
prm.declare_entry("lambda austenite", "144.0e9",
Patterns::Double(),
"lambda austenite");
prm.declare_entry("mu austenite", "74.0e9",
Patterns::Double(0.0),
"mu austenite");
prm.declare_entry("lambda martensite", "379.0e9",
Patterns::Double(),
"lambda martensite");
prm.declare_entry("mu martensite", "134.0e9",
Patterns::Double(0.0),
"mu martensite");
prm.declare_entry("kinetic coeff", "2600.0",
Patterns::Double(0.0),
"kinetic coeff");
prm.declare_entry("energy coeff", "2.59e-10",
Patterns::Double(0.0),
"energy coeff");
prm.declare_entry("barrier height", "28.0e6",
Patterns::Double(),
"barrier height");
prm.declare_entry("temperature", "50.0",
Patterns::Double(0.0),
"temperature");
prm.declare_entry("temperature crit", "-183.0",
Patterns::Double(),
"temperature crit");
prm.declare_entry("temperature eql", "215.0",
Patterns::Double(),
"temperature eql");
prm.declare_entry("crank-nicolson parameter", "0.5",
Patterns::Double(),
"crank-nicolson parameter");
}
prm.leave_subsection();
}
void Materials::parse_parameters(ParameterHandler &prm)
{
prm.enter_subsection("Material properties");
{
lambda0 = prm.get_double("lambda austenite");
mu0 = prm.get_double("mu austenite");
lambda1 = prm.get_double("lambda martensite");
mu1 = prm.get_double("mu martensite");
L = prm.get_double("kinetic coeff");
beta = prm.get_double("energy coeff");
A0 = prm.get_double("barrier height");
theta = prm.get_double("temperature");
thetac = prm.get_double("temperature crit");
thetae = prm.get_double("temperature eql");
thetan = prm.get_double("crank-nicolson parameter");
}
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", "2e-12",
Patterns::Double(),
"End time");
prm.declare_entry("Time step size", "1e-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 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);
Time::declare_parameters(prm);
}
void AllParameters::parse_parameters(ParameterHandler &prm)
{
FESystem::parse_parameters(prm);
Geometry::parse_parameters(prm);
Materials::parse_parameters(prm);
Time::parse_parameters(prm);
}
}
// DEFINE SECOND ORDER IDENTITY, AND TWO FOURTH ORDER IDENTITY TENSORS
template <int dim>
class StandardTensors
{
public:
static const SymmetricTensor<2, dim> I;
static const SymmetricTensor<4, dim> IxI;
static const SymmetricTensor<4, dim> II;
// static const SymmetricTensor<2, dim> transformation_strain;
};
template <int dim>
const SymmetricTensor<2, dim>
StandardTensors<dim>::I = unit_symmetric_tensor<dim>();
template <int dim>
const SymmetricTensor<4, dim>
StandardTensors<dim>::IxI = outer_product(I, I);
template <int dim>
const SymmetricTensor<4, dim>
StandardTensors<dim>::II = identity_tensor<dim>();
// DEFINE TIME STEP, CURRENT TIME ETC.
class Time
{
public:
Time (const double time_end,
const double delta_t)
:
timestep(0),
time_current(0.0),
time_end(time_end),
delta_t(delta_t)
{}
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_Compressible_Neo_Hook_Three_Field
{
public:
Material_Compressible_Neo_Hook_Three_Field(const double lambda0,
const double mu0,
const double lambda1,
const double mu1,
const double A0,
const double theta,
const double thetac,
const double thetae)
:
det_F(1.0),
I1(0.0), // I1 = trace(Be)
ge(StandardTensors<dim>::I), // = Fe.Fe^T
Be(SymmetricTensor<2, dim>()), // = 0.5(Fe.Fe^T-I)
lambda00(lambda0),
mu00(mu0),
lambda11(lambda1),
mu11(mu1),
A00(A0),
theta1(theta),
thetac1(thetac),
thetae1(thetae)
{}
~Material_Compressible_Neo_Hook_Three_Field()
{}
void update_material_data (const Tensor<2, dim> &F, const Tensor<2, dim> &Fe,
const double phi, const double dphi,const double ddphi,
const double dPotentialBarrier, const double ddPotentialBarrier)
{
det_F = determinant(F); // determinant of total F
ge = symmetrize(Fe*transpose(Fe)); // ge = Fe.Fe^T
Be = 0.5*(ge-StandardTensors<dim>::I); // Be = 1/2(ge-I)
I1 = trace(Be); // first invariant of Be
lambda = lambda00+(lambda11-lambda00)*phi; // definition of the Lame constants are obvious
mu = mu00+(mu11-mu00)*phi;
deltaG = A00*(theta1-thetae1)/3.0; // difference in thermal energy
A = A00*(theta1-thetac1); // compute barrier height
dphi1 = dphi; // d_phi/d_eta
ddphi1 = ddphi;
dPotentialBarrier1 = dPotentialBarrier;
ddPotentialBarrier1 = ddPotentialBarrier;
Assert(det_F > 0, ExcInternalError());
}
SymmetricTensor<2, dim> get_tau() // tau is the Kirchhoff stress := J*Cauchy stress
{
return get_tau_kirchhoff();
}
SymmetricTensor<4, dim> get_Jc() const // comupte the fourth order modulus tensor in the deformed config.
{
return get_Jc_modulus();
}
double get_det_F() const
{
return det_F;
}
double get_driving_force_noStress() const // driving force excluding the transformational work
{
return get_driv_forc_noStress ();
}
double get_deri_driving_force_noStress() const // driving force excluding the transformational work
{
return get_deri_driv_forc_noStress ();
}
protected:
double det_F;
double I1;
SymmetricTensor<2, dim> ge;
SymmetricTensor<2, dim> Be;
double lambda00;
double mu00;
double lambda11;
double mu11;
double A00;
double theta1;
double thetac1;
double thetae1;
double lambda;
double mu;
double deltaG;
double A;
double dphi1;
double ddphi1;
double dPotentialBarrier1;
double ddPotentialBarrier1;
// compute the Kirchhoff stress
SymmetricTensor<2, dim> get_tau_kirchhoff() const
{
SymmetricTensor<2, dim> kirchhoff_stress = lambda*I1*ge +
mu*symmetrize(Tensor<2, dim>(ge) * Tensor<2, dim>(ge))-mu*ge;
return kirchhoff_stress;
}
// compute the modulus J*d_ijkl
SymmetricTensor<4, dim> get_Jc_modulus() const
{
SymmetricTensor<4,dim> elasticityTensor;
for (unsigned int i=0; i<dim; ++i)
for (unsigned int j=0; j<dim; ++j)
for (unsigned int k=0; k<dim; ++k)
for (unsigned int l=0; l<dim; ++l)
elasticityTensor[i][j][k][l] =
lambda*ge[i][j]*ge[k][l]+mu*(ge[i][l]*ge[j][k]+ge[i][k]*ge[j][l]);
return elasticityTensor;
}
// compute the driving force excluding the transformational work
double get_driv_forc_noStress () const
{
return -A*dPotentialBarrier1-deltaG*dphi1;
}
double get_deri_driv_forc_noStress () const
{
return -A*ddPotentialBarrier1-deltaG*ddphi1;
}
};
// updates the quadrature point history
template <int dim>
class PointHistory
{
public:
PointHistory()
:
material(NULL),
F_inv(StandardTensors<dim>::I),
tau(SymmetricTensor<2, dim>()),
Jc(SymmetricTensor<4, dim>())
{}
virtual ~PointHistory()
{
delete material;
material = NULL;
}
void setup_lqp (const Parameters::AllParameters ¶meters)
{
material = new Material_Compressible_Neo_Hook_Three_Field<dim>(parameters.lambda0,
parameters.mu0, parameters.lambda1, parameters.mu1, parameters.A0,
parameters.theta, parameters.thetac, parameters.thetae);
update_values(Tensor<2, dim>(), double (), double (), double());
}
void update_values (const Tensor<2, dim> &Grad_u_n,
const double new_eta,
const double old_eta,
const double thetan)
{
double eta = thetan*new_eta +(1-thetan)*old_eta;
const double phi = 3*eta*eta-2*eta*eta*eta; // interpolation function
const double dphi = 6*eta-6*eta*eta; // derivative of interpolation fn
const double ddphi = 6-12*eta;
const double dPotentialBarrier = 2*eta-6*eta*eta+4*eta*eta*eta;
const double ddPotentialBarrier =2-12*eta+12*eta*eta;
// double strectch1 = 0.8;
// double strectch2 = 1.2;
Tensor<2, dim> tmp; // transformation strain tensor
// tmp[0][0] = pow(strectch1, 1.5)/(1.0+strectch1) + pow(strectch2, 1.5)/(1.0+strectch2)-1.0;
// tmp[1][1] = pow(strectch1, 0.5)/(1.0+strectch1) + pow(strectch2, 0.5)/(1.0+strectch2)-1.0;
// tmp[2][2] = 0.0;
// tmp[0][1] = -strectch1/(1.0+strectch1) + strectch2/(1.0+strectch2);
// tmp[1][0] = -strectch1/(1.0+strectch1) + strectch2/(1.0+strectch2);
// tmp[0][2] = 0.0;
// tmp[2][0] = 0.0;
// tmp[1][2] = 0.0;
// tmp[2][1] = 0.0;
//
tmp[0][0] = 0.1;
tmp[1][1] = -0.05;
tmp[2][2] = -0.05;
tmp[0][1] = 0.0;
tmp[1][0] = 0.0;
tmp[0][2] = 0.0;
tmp[2][0] = 0.0;
tmp[1][2] = 0.0;
tmp[2][1] = 0.0;
SymmetricTensor<2, dim> transStrain = symmetrize(tmp);
const Tensor<2, dim> Ft = (Tensor<2, dim>(StandardTensors<dim>::I) +
Tensor<2, dim>(transStrain)*phi); // transformation deformation gradient
F = (Tensor<2, dim>(StandardTensors<dim>::I) + Grad_u_n); // total F
Fe = F * invert(Ft); // elastic part of deformation grad
material->update_material_data(F, Fe, phi, dphi, ddphi, dPotentialBarrier, ddPotentialBarrier);
F_inv = invert(F);
tau = material->get_tau(); // extracting kirchhoff stress
Jc = material->get_Jc(); // extracting J*d_ijkl
driving_force_noStress = material->get_driving_force_noStress(); // extracting driving force with no stress
deri_driving_force_noStress = material->get_deri_driving_force_noStress();
const Tensor<2, dim> temp_tensor = F_inv*Tensor<2, dim>(tau);
const Tensor<2, dim> temp_tensor1 = temp_tensor * Fe;
get_driv_forc = scalar_product(temp_tensor1, tmp)*dphi + driving_force_noStress; // computing total driving force
get_deri_driv_forc = scalar_product(temp_tensor1, tmp)*ddphi + deri_driving_force_noStress;
Assert(determinant(F_inv) > 0, ExcInternalError());
}
const Tensor<2, dim> &get_F() const
{
return F;
}
double get_det_F() const
{
return material->get_det_F();
}
const Tensor<2, dim> &get_F_inv() const
{
return F_inv;
}
const SymmetricTensor<2, dim> &get_tau() const
{
return tau;
}
const SymmetricTensor<4, dim> &get_Jc() const
{
return Jc;
}
double get_driving_force() const
{
return get_driv_forc;
}
double get_deri_driving_force() const
{
return get_deri_driv_forc;
}
private:
Material_Compressible_Neo_Hook_Three_Field<dim> *material;
Tensor<2, dim> F_inv;
Tensor<2, dim> F;
SymmetricTensor<2, dim> tau;
SymmetricTensor<4, dim> Jc;
Tensor<2, dim> Fe;
double driving_force_noStress;
double deri_driving_force_noStress;
double dphi;
double ddphi;
double get_driv_forc;
double get_deri_driv_forc;
};
///////////////////////////////////////////////////////////////
template <int dim>
class Solid
{
public:
Solid(const std::string &input_file);
virtual
~Solid();
void
run();
private:
void make_grid();
void system_setup();
void assemble_system();
void make_constraints(const int &it_nr);
void solve_nonlinear_timestep();
unsigned int solve();
void assemble_system_eta();
void solve_nonlinear_timestep_eta();
unsigned int solve_eta();
void setup_qph();
void update_qph_incremental();
void output_results() const;
Parameters::AllParameters parameters;
double vol_reference; // volume of the reference config
double vol_current; // volume of the current config
Time time; // variable of type class 'Time'
TimerOutput timer;
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<dim> triangulation;
ConditionalOStream pcout;
const unsigned int degree; // degree of polynomial of shape functions
const FESystem<dim> fe; // fe object
DoFHandler<dim> dof_handler; // we have used two dof_handler: one for mechanics another for order parameter
const unsigned int dofs_per_cell; // no of dofs per cell for the mechanics problem
const FEValuesExtractors::Vector u_fe;
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
FE_DGQ<dim> history_fe;
DoFHandler<dim> history_dof_handler;
std::vector<PointHistory<dim> > quadrature_point_history;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
LA::MPI::SparseMatrix tangent_matrix; // tangent stiffenss matrix
LA::MPI::Vector system_rhs; // system right hand side or residual of mechanics problem
LA::MPI::Vector solution; // solution vector for displacement
LA::MPI::Vector solution_update; // another vector containing the displacement soln
LA::MPI::Vector tmp;
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 QGauss<dim> qf_cell_eta;
const unsigned int n_q_points_eta;
IndexSet locally_owned_dofs_eta;
IndexSet locally_relevant_dofs_eta;
LA::MPI::SparseMatrix mass_matrix; // mass matrix out of Ginxburg-Landau eqn
LA::MPI::SparseMatrix laplace_matrix; // Laplace matrix out of Ginxburg-Landau eqn
LA::MPI::SparseMatrix system_matrix_eta;
LA::MPI::SparseMatrix nl_matrix;
LA::MPI::SparseMatrix tmp_matrix;
LA::MPI::Vector nl_term;
LA::MPI::Vector system_rhs_eta;
LA::MPI::Vector solution_eta; // solution vector for eta
LA::MPI::Vector old_solution_eta;
LA::MPI::Vector solution_update_eta;
LA::MPI::Vector tmp_vector; // a vector used for solving eta
LA::MPI::Vector tmp_eta;
};
/////////////////////////////////////////////////////////
// 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
{
// return 0.01;
return (rand()%100+0.0)*0.01;
// 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 std::max(tanh((1e-9-sqrt((p[0]-5.0e-9)*(p[0]-5.0e-9)+(p[1]-5.0e-9)*(p[1]-5.0e-9)))/0.4e-9),0.0);
//return (0.5*std::tanh(3.0*(p[0]-5.0e-9)/7.07e-10)+0.5);
//return 0.1 * std::max(tanh((p[0]-5.0e-9)/(0.4e-9)), 0.0);
// return 0.5*std::tanh(p[0]/(0.4e-9))+0.5;
}
/////////////////////////////////////////////////
template <int dim>
class BoundaryDisplacement : public Function<dim>
{
public:
BoundaryDisplacement (const int direction, const int time_step);
virtual
void
vector_value (const Point<dim> &p,
Vector<double> &values) const;
virtual
void
vector_value_list (const std::vector<Point<dim> > &points,
std::vector<Vector<double> > &value_list) const;
private:
const int direction;
const int time_step;
};
template <int dim>
BoundaryDisplacement<dim>::BoundaryDisplacement (const int direction,const int time_step)
:
Function<dim> (dim),
direction(direction),
time_step(time_step)
{}
template <int dim>
inline
void
BoundaryDisplacement<dim>::vector_value (const Point<dim> &/*p*/,
Vector<double> &values) const
{
Assert (values.size() == dim,
ExcDimensionMismatch (values.size(), dim));
values = 0;
if (time_step <3)
values(direction) = 0.1e-9;
else
values(direction)= 0.001e-9;
}
template <int dim>
void
BoundaryDisplacement<dim>::vector_value_list (const std::vector<Point<dim> > &points,
std::vector<Vector<double> > &value_list) const
{
const unsigned int n_ooints = points.size();
Assert (value_list.size() == n_ooints,
ExcDimensionMismatch (value_list.size(), n_ooints));
for (unsigned int p=0; p<n_ooints; ++p)
BoundaryDisplacement<dim>::vector_value (points[p],
value_list[p]);
}
///////////////////////////////////////////////////
// constructor
template <int dim>
Solid<dim>::Solid(const std::string &input_file)
:
mpi_communicator (MPI_COMM_WORLD),
triangulation (mpi_communicator,
typename Triangulation<dim>::MeshSmoothing
(Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening)),
parameters(input_file),
time(parameters.end_time, parameters.delta_t),
pcout (std::cout,
(Utilities::MPI::this_mpi_process(mpi_communicator)
== 0)),
timer(mpi_communicator,
pcout,
TimerOutput::summary,
TimerOutput::wall_times),
degree(parameters.poly_degree),
fe(FE_Q<dim>(parameters.poly_degree), dim), // displacement
dof_handler(triangulation),
dofs_per_cell (fe.dofs_per_cell),
u_fe(0),
qf_cell(parameters.quad_order),
qf_face(parameters.quad_order),
n_q_points (qf_cell.size()),
n_q_points_f (qf_face.size()),
degree_eta(parameters.poly_degree),
fe_eta (parameters.poly_degree),
dof_handler_eta (triangulation),
dofs_per_cell_eta(fe_eta.dofs_per_cell),
qf_cell_eta(parameters.quad_order),
n_q_points_eta (qf_cell_eta.size()),
history_dof_handler (triangulation),
history_fe (parameters.poly_degree)
{}
//destructor
template <int dim>
Solid<dim>::~Solid()
{
dof_handler.clear();
dof_handler_eta.clear();
}
//////////////////////////////////////////////////
template <int dim>
void Solid<dim>::run()
{
make_grid(); // generates the geometry and mesh
system_setup(); // sets up the system matrices
VectorTools::interpolate(dof_handler_eta, InitialValues<dim>(), tmp_eta); //initial eta
solution_eta= tmp_eta;
old_solution_eta = solution_eta;
update_qph_incremental();
// this is the zeroth iteration for compute the initial distorted solution
//of the body due to arbitrary distribution of initial eta
solve_nonlinear_timestep();
output_results();
time.increment();
// computed actual time integration to update displacement and eta
while (time.current() <= time.end())
{
pcout << std::endl
<< "Time step #" << time.get_timestep() << "; "
<< "advancing to t = " << time.current() << "."
<< std::endl;
solve_nonlinear_timestep_eta();
solve_nonlinear_timestep();
output_results();
time.increment();
}
}
///////////////////////////////////
template <int dim>
void Solid<dim>::solve_nonlinear_timestep()
{
double initial_rhs_norm = 0.;
unsigned int newton_iteration = 0;
for (; newton_iteration < 20; ++newton_iteration)
{
assemble_system ();
make_constraints(newton_iteration);
if (newton_iteration == 0)
initial_rhs_norm = system_rhs.l2_norm();
const unsigned int n_iterations
= solve ();
tmp += solution_update;
solution = tmp;
update_qph_incremental();
if (newton_iteration == 0)
std::cout << " Solving for Displacement: " << n_iterations;
else
std::cout << '+' << n_iterations;
if (newton_iteration > 0 && system_rhs.l2_norm() <= 1e-6 * initial_rhs_norm)
{
std::cout << " CG iterations per nonlinear step. CONVERGED! " << std::endl;
break;
}
AssertThrow (newton_iteration < 19,
ExcMessage("No convergence in nonlinear solver!"));
}
}
/////////////////////////////////////////
template <int dim>
void Solid<dim>::solve_nonlinear_timestep_eta()
{
old_solution_eta = solution_eta;
double initial_rhs_norm = 0.;
unsigned int newton_iteration = 0;
for (; newton_iteration < 20; ++newton_iteration)
{
assemble_system_eta ();
if (newton_iteration == 0)
initial_rhs_norm = system_rhs_eta.l2_norm();
const unsigned int n_iterations
= solve_eta ();
tmp_eta += solution_update_eta;
solution_eta= tmp_eta;
update_qph_incremental();
if (newton_iteration == 0)
std::cout << " Solving for Order Parameter: " << n_iterations;
else
std::cout << '+' << n_iterations;
if (newton_iteration > 0 && system_rhs_eta.l2_norm() <= 1e-6 * initial_rhs_norm)
{
std::cout << " CG iterations per nonlinear step. CONVERGED! " << std::endl;
break;
}
AssertThrow (newton_iteration < 19,
ExcMessage("No convergence in nonlinear solver!"));
}
}
///////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::make_grid()
{
// // Divide the beam, but only along the x- and y-coordinate directions
// std::vector< unsigned int > repetitions(dim, 32);
// // Only allow one element through the thickness
// // (modelling a plane strain condition)
// if (dim == 3)
// repetitions[dim-1] = 1;
//
// const Point<dim> bottom_left = (dim == 3 ? Point<dim>(0.0, 0.0, -0.5) : Point<dim>(0.0, 0.0));
// const Point<dim> top_right = (dim == 3 ? Point<dim>(10.0, 10.0, 0.5) : Point<dim>(10.0, 10.0));
//
// GridGenerator::subdivided_hyper_rectangle(triangulation,
// repetitions,
// bottom_left,
// top_right,
// true);
GridGenerator::hyper_rectangle(triangulation,
Point<dim>(0.0, 0.0, 0.0),
Point<dim>(10.0, 10.0, 10.0),
true);
triangulation.refine_global (parameters.refinement);
GridTools::scale(parameters.scale, triangulation);
}
///////////////////////////////////////////////////
template <int dim>
void Solid<dim>::system_setup()
{
dof_handler.distribute_dofs(fe);
dof_handler_eta.distribute_dofs (fe_eta);
history_dof_handler.distribute_dofs (history_fe);
const unsigned int n_dofs = dof_handler.n_dofs(),
n_dofs_eta = dof_handler_eta.n_dofs();
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Total number of cells: "
<< triangulation.n_cells()
<< std::endl
<< " Number of degrees of freedom: "
<< n_dofs + n_dofs_eta
<< " (" << n_dofs << '+' << n_dofs_eta << ')'
<< std::endl;
locally_owned_dofs = dof_handler.locally_owned_dofs ();
DoFTools::extract_locally_relevant_dofs (dof_handler , locally_relevant_dofs);
DynamicSparsityPattern dsp(locally_relevant_dofs);
DoFTools::make_sparsity_pattern (dof_handler, dsp);
SparsityTools::distribute_sparsity_pattern (dsp,
dof_handler.n_locally_owned_dofs_per_processor(),
mpi_communicator,
locally_relevant_dofs);
tangent_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
solution_update.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
system_rhs.reinit(locally_owned_dofs,
mpi_communicator);
tmp.reinit(locally_owned_dofs,
mpi_communicator);
locally_owned_dofs_eta = dof_handler_eta.locally_owned_dofs ();
DoFTools::extract_locally_relevant_dofs (dof_handler_eta , locally_relevant_dofs_eta);
DynamicSparsityPattern dsp_eta(locally_relevant_dofs_eta);
DoFTools::make_sparsity_pattern (dof_handler_eta, dsp_eta);
SparsityTools::distribute_sparsity_pattern (dsp_eta,
dof_handler_eta.n_locally_owned_dofs_per_processor(),
mpi_communicator,
locally_relevant_dofs_eta);
mass_matrix.reinit (locally_owned_dofs_eta,
locally_owned_dofs_eta,
dsp_eta,
mpi_communicator);
laplace_matrix.reinit (locally_owned_dofs_eta,
locally_owned_dofs_eta,
dsp_eta,
mpi_communicator);
system_matrix_eta.reinit (locally_owned_dofs_eta,
locally_owned_dofs_eta,
dsp_eta,
mpi_communicator);
nl_matrix.reinit (locally_owned_dofs_eta,
locally_owned_dofs_eta,
dsp_eta,
mpi_communicator);
tmp_matrix.reinit (locally_owned_dofs_eta,
locally_owned_dofs_eta,
dsp_eta,
mpi_communicator);
// MatrixCreator::create_mass_matrix(dof_handler_eta,
// QGauss<dim>(fe_eta.degree+1),
// mass_matrix);
// MatrixCreator::create_laplace_matrix(dof_handler_eta,
// QGauss<dim>(fe_eta.degree+1),
// laplace_matrix);
solution_eta.reinit(locally_owned_dofs_eta,
locally_relevant_dofs_eta,
mpi_communicator);
old_solution_eta.reinit(locally_owned_dofs_eta,
locally_relevant_dofs_eta,
mpi_communicator);
solution_update_eta.reinit(locally_owned_dofs_eta,
locally_relevant_dofs_eta,
mpi_communicator);
system_rhs_eta.reinit(locally_owned_dofs_eta,
mpi_communicator);
nl_term.reinit(locally_owned_dofs_eta,
mpi_communicator);
tmp_vector.reinit(locally_owned_dofs_eta,
mpi_communicator);
tmp_eta.reinit(locally_owned_dofs_eta,
mpi_communicator);
// sets up the quadrature point history
setup_qph();
}
//////////////////////////////////
template <int dim>
void Solid<dim>::make_constraints(const int &it_nr)
{
if (it_nr > 1)
return;
constraints.clear();
constraints.reinit (locally_relevant_dofs);
const bool apply_dirichlet_bc = (it_nr == 0);
const int timestep = time.get_timestep();
const FEValuesExtractors::Scalar x_displacement(0);
const FEValuesExtractors::Scalar y_displacement(1);
const FEValuesExtractors::Scalar z_displacement(2);
{
const int boundary_id = 0;
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement)/*|
fe.component_mask(y_displacement) |
fe.component_mask(z_displacement)*/);
else
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement)/* |
fe.component_mask(y_displacement) |
fe.component_mask(z_displacement)*/);
}
{
const int boundary_id = 1;
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
BoundaryDisplacement<dim>(0, timestep),
constraints,
fe.component_mask(x_displacement));
else
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
}
//
////
{
const int boundary_id = 2;
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(y_displacement));
else
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(y_displacement));
}
//
// {
// const int boundary_id = 3;
//
// if (apply_dirichlet_bc == true)
// VectorTools::interpolate_boundary_values(dof_handler,
// boundary_id,
// ZeroFunction<dim>(dim),
// constraints,
// fe.component_mask(y_displacement));
// else
// VectorTools::interpolate_boundary_values(dof_handler,
// boundary_id,
// ZeroFunction<dim>(dim),
// constraints,
// fe.component_mask(y_displacement));
// }
{
const int boundary_id = 4;
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(z_displacement));
else
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(z_displacement));
}
// {
// const int boundary_id = 5;
//
// if (apply_dirichlet_bc == true)
// VectorTools::interpolate_boundary_values(dof_handler,
// boundary_id,
// ZeroFunction<dim>(dim),
// constraints,
// fe.component_mask(z_displacement));
// else
// VectorTools::interpolate_boundary_values(dof_handler,
// boundary_id,
// ZeroFunction<dim>(dim),
// constraints,
// fe.component_mask(z_displacement));
// }
constraints.close();
}
/////////////////////////////////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::assemble_system ()
{
tangent_matrix = 0;
system_rhs = 0;
FEValues<dim> fe_values (fe, qf_cell,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEFaceValues<dim> fe_face_values (fe, qf_face,
update_values | update_quadrature_points |
update_normal_vectors | update_JxW_values);
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
std::vector<double> Nx(dofs_per_cell);
std::vector<Tensor<2, dim> > grad_Nx(dofs_per_cell);
std::vector<SymmetricTensor<2, dim> > symm_grad_Nx(dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
{
fe_values.reinit (cell);
cell_matrix = 0;
cell_rhs = 0;
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 Tensor<2, dim> tau = lqph[q_point].get_tau();
const SymmetricTensor<2, dim> symm_tau = lqph[q_point].get_tau();
const SymmetricTensor<4, dim> Jc = lqph[q_point].get_Jc();
const double JxW = fe_values.JxW(q_point);
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
grad_Nx[k] = fe_values[u_fe].gradient(k, q_point) * F_inv;
symm_grad_Nx[k] = symmetrize(grad_Nx[k]);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int component_i = fe.system_to_component_index(i).first;
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const unsigned int component_j = fe.system_to_component_index(j).first;
cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution:
* symm_grad_Nx[j] * JxW;
if (component_i == component_j) // geometrical stress contribution
cell_matrix(i, j) += grad_Nx[i][component_i] * tau
* grad_Nx[j][component_j] * JxW;
}
cell_rhs(i) -= symm_grad_Nx[i] * symm_tau * JxW;
}
}
//
//
// for (unsigned int face_number=0; face_number<GeometryInfo<dim>::faces_per_cell; ++face_number)
// if (cell->face(face_number)->at_boundary()
// &&
// (cell->face(face_number)->boundary_id() == 1))
// {
// fe_face_values.reinit (cell, face_number);
// for (unsigned int q_point=0; q_point<n_q_points_f; ++q_point)
// {
// double magnitude;
//
// if (cell->face(face_number)->center()[1]>= -0.2e-9 && cell->face(face_number)->center()[1]<= 0.2e-9)
// magnitude = -10e8;
// else
// magnitude = 0.0;
//
// static const Tensor<1, dim> dir ({1.0,0.0,0.0});
// const Tensor<1, dim> traction = magnitude*dir;
//
//
// for (unsigned int i=0; i<dofs_per_cell; ++i)
// {
// const unsigned int
// component_i = fe.system_to_component_index(i).first;
// cell_rhs(i) += (traction[component_i] *
// fe_face_values.shape_value(i,q_point) *
// fe_face_values.JxW(q_point));
// }
// }
// }
//
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
tangent_matrix,
system_rhs);
}
tangent_matrix.compress (VectorOperation::add);
system_rhs.compress (VectorOperation::add);
}
//////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::assemble_system_eta ()
{
system_matrix_eta = 0;
system_rhs_eta = 0;
mass_matrix = 0;
laplace_matrix = 0;
nl_matrix = 0;
nl_term = 0;
FEValues<dim> fe_values_eta (fe_eta, qf_cell_eta,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FullMatrix<double> cell_mass_matrix (dofs_per_cell_eta, dofs_per_cell_eta);
FullMatrix<double> cell_laplace_matrix (dofs_per_cell_eta, dofs_per_cell_eta);
FullMatrix<double> cell_nl_matrix (dofs_per_cell_eta, dofs_per_cell_eta);
Vector <double> cell_nl_term (dofs_per_cell_eta);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell_eta);
std::vector<double> phi (dofs_per_cell_eta);
std::vector<Tensor<1,dim> > grad_phi (dofs_per_cell_eta);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler_eta.begin_active(),
endc = dof_handler_eta.end();
for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
{
fe_values_eta.reinit(cell);
cell_mass_matrix = 0;
cell_laplace_matrix = 0;
cell_nl_matrix = 0;
cell_nl_term = 0;
PointHistory<dim> *lqph =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
for (unsigned int q=0; q<n_q_points_eta; ++q)
{
const double driving_force = lqph[q].get_driving_force();
const double deri_driving_force = lqph[q].get_deri_driving_force();
for (unsigned int k=0; k<dofs_per_cell_eta; ++k)
{
phi[k] = fe_values_eta.shape_value (k, q);
grad_phi[k] = fe_values_eta.shape_grad (k, q);
}
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) += (phi[i] * phi[j]) * fe_values_eta.JxW(q);
cell_laplace_matrix(i,j) += (grad_phi[i] * grad_phi[j]) * fe_values_eta.JxW(q);
cell_nl_matrix(i,j) += deri_driving_force * phi[i] * phi[j] * fe_values_eta.JxW(q);
}
cell_nl_term(i) += driving_force * phi[i] * fe_values_eta.JxW (q);
}
}
cell->get_dof_indices (local_dof_indices);
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[i],
local_dof_indices[j],
cell_mass_matrix(i,j));
laplace_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_laplace_matrix(i,j));
nl_matrix. add (local_dof_indices[i],
local_dof_indices[j],
cell_nl_matrix(i,j));
}
nl_term (local_dof_indices[i]) += cell_nl_term(i);
}
}
mass_matrix.compress (VectorOperation::add);
laplace_matrix.compress (VectorOperation::add);
nl_matrix.compress (VectorOperation::add);
nl_term.compress (VectorOperation::add);
system_matrix_eta.copy_from (mass_matrix);
system_matrix_eta.add ( laplace_matrix, parameters.L * parameters.beta * parameters.delta_t * parameters.thetan);
system_matrix_eta.add ( nl_matrix, -parameters.L * parameters.delta_t * parameters.thetan);
tmp_matrix.copy_from (mass_matrix);
tmp_matrix.add ( laplace_matrix, parameters.L * parameters.beta * parameters.delta_t * parameters.thetan);
tmp_matrix.vmult (tmp_vector, solution_eta);
system_rhs_eta += tmp_vector;
tmp_matrix.add( laplace_matrix, -parameters.L * parameters.beta * parameters.delta_t);
tmp_matrix.vmult (tmp_vector, old_solution_eta);
system_rhs_eta -= tmp_vector;
system_rhs_eta.add (-parameters.L * parameters.delta_t, nl_term);
system_rhs_eta *= -1;
system_matrix_eta.compress (VectorOperation::add);
system_rhs_eta.compress (VectorOperation::add);
}
///////////////////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::setup_qph()
{
{
unsigned int our_cells = 0;
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end(); ++cell)
if (cell->is_locally_owned())
++our_cells;
triangulation.clear_user_data();
{
std::vector<PointHistory<dim> > tmp;
tmp.swap (quadrature_point_history);
}
quadrature_point_history.resize (our_cells * n_q_points);
unsigned int history_index = 0;
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end(); ++cell)
if (cell->is_locally_owned())
{
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)
if (cell->is_locally_owned())
{
PointHistory<dim> *lqph =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
Assert(lqph >= &quadrature_point_history.front(), ExcInternalError());
Assert(lqph <= &quadrature_point_history.back(), ExcInternalError());
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
lqph[q_point].setup_lqp(parameters);
}
}
// Updates the data at all quadrature points over a loop run by WorkStream::run
template <int dim>
void Solid<dim>::update_qph_incremental()
{
FEValues<dim> fe_values (fe, qf_cell,
update_values | update_gradients);
FEValues<dim> fe_values_eta (fe_eta, qf_cell,
update_values | update_gradients);
std::vector<Tensor<2, dim> > solution_grads_total (qf_cell.size());
std::vector<double> solution_eta_total (qf_cell.size());
std::vector<double> old_solution_eta_total (qf_cell.size());
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
typename DoFHandler<dim>::active_cell_iterator
cell_eta = dof_handler_eta.begin_active();
for (; cell!=endc; ++cell, ++cell_eta)
if (cell->is_locally_owned())
{
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(solution_grads_total.size() == n_q_points,
ExcInternalError());
Assert(solution_eta_total.size() == n_q_points,
ExcInternalError());
fe_values.reinit(cell);
fe_values_eta.reinit (cell_eta);
fe_values[u_fe].get_function_gradients(solution, solution_grads_total);
fe_values_eta.get_function_values(solution_eta, solution_eta_total);
fe_values_eta.get_function_values(old_solution_eta, old_solution_eta_total);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
lqph[q_point].update_values(solution_grads_total[q_point],
solution_eta_total[q_point],
old_solution_eta_total[q_point],
parameters.thetan);
}
}
//////////////////////////////
template <int dim>
unsigned int
Solid<dim>::solve ()
{
LA::MPI::Vector
completely_distributed_solution (locally_owned_dofs, mpi_communicator);
SolverControl solver_control (dof_handler.n_dofs(), 1e-12*system_rhs.l2_norm());
#ifdef USE_PETSC_LA
LA::SolverCG solver(solver_control, mpi_communicator);
#else
LA::SolverCG solver(solver_control);
#endif
LA::MPI::PreconditionAMG preconditioner;
LA::MPI::PreconditionAMG::AdditionalData data;
#ifdef USE_PETSC_LA
data.symmetric_operator = true;
//data.strong_threshold = 0.5;
#else
/* Trilinos defaults are good */
#endif
preconditioner.initialize(tangent_matrix, data);
// LA::MPI::PreconditionJacobi preconditioner(tangent_matrix);
solver.solve (tangent_matrix, completely_distributed_solution, system_rhs,
preconditioner);
constraints.distribute (completely_distributed_solution);
solution_update = completely_distributed_solution;
return solver_control.last_step();
}
///////////////////////////////////////////////////////////////////////////
template <int dim>
unsigned int
Solid<dim>::solve_eta ()
{
LA::MPI::Vector
completely_distributed_solution_eta (locally_owned_dofs_eta, mpi_communicator);
SolverControl solver_control (dof_handler_eta.n_dofs(), 1e-12*system_rhs_eta.l2_norm());
#ifdef USE_PETSC_LA
LA::SolverCG solver(solver_control, mpi_communicator);
#else
LA::SolverCG solver(solver_control);
#endif
LA::MPI::PreconditionAMG preconditioner;
LA::MPI::PreconditionAMG::AdditionalData data;
#ifdef USE_PETSC_LA
data.symmetric_operator = true;
#else
/* Trilinos defaults are good */
#endif
preconditioner.initialize(system_matrix_eta, data);
//LA::MPI::PreconditionJacobi preconditioner(system_matrix_eta);
solver.solve (system_matrix_eta, completely_distributed_solution_eta, system_rhs_eta,
preconditioner);
solution_update_eta = completely_distributed_solution_eta;
return solver_control.last_step();
}
///////////////////////////////////////////////////////////
template <int dim>
void Solid<dim>::output_results() const
{
DataOut<dim> data_out;
std::vector<std::string> displacement_names;
switch (dim)
{
case 1:
displacement_names.push_back ("displacement");
break;
case 2:
displacement_names.push_back ("x_displacement");
displacement_names.push_back ("y_displacement");
break;
case 3:
displacement_names.push_back ("x_displacement");
displacement_names.push_back ("y_displacement");
displacement_names.push_back ("z_displacement");
break;
default:
Assert (false, ExcNotImplemented());
}
data_out.add_data_vector (dof_handler, solution, displacement_names);
data_out.add_data_vector (dof_handler_eta, solution_eta, "order_parameter");
/////////////////////////
Vector<double> norm_of_stress (triangulation.n_active_cells());
{
typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
{
SymmetricTensor<2,dim> accumulated_stress;
for (unsigned int q=0; q<qf_cell.size(); ++q)
accumulated_stress +=
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer())[q].get_tau();
norm_of_stress(cell->active_cell_index())
= (accumulated_stress /
qf_cell.size()).norm();
}
else
norm_of_stress(cell->active_cell_index()) = -1e+20;
}
data_out.add_data_vector (norm_of_stress, "norm_of_stress");
std::vector< std::vector< Vector<double> > >
history_field (dim, std::vector< Vector<double> >(dim)),
local_history_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
local_history_fe_values (dim, std::vector< Vector<double> >(dim)),
history_field_stress (dim, std::vector< Vector<double> >(dim)),
local_history_values_at_qpoints_stress (dim, std::vector< Vector<double> >(dim)),
local_history_fe_values_stress (dim, std::vector< Vector<double> >(dim));
for (unsigned int i=0; i<dim; i++)
for (unsigned int j=0; j<dim; j++)
{
history_field[i][j].reinit(history_dof_handler.n_dofs());
local_history_values_at_qpoints[i][j].reinit(qf_cell.size());
local_history_fe_values[i][j].reinit(history_fe.dofs_per_cell);
history_field_stress[i][j].reinit(history_dof_handler.n_dofs());
local_history_values_at_qpoints_stress[i][j].reinit(qf_cell.size());
local_history_fe_values_stress[i][j].reinit(history_fe.dofs_per_cell);
}
FullMatrix<double> qpoint_to_dof_matrix (history_fe.dofs_per_cell,
qf_cell.size());
FETools::compute_projection_from_quadrature_points_matrix
(history_fe,
qf_cell, qf_cell,
qpoint_to_dof_matrix);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end(),
dg_cell = history_dof_handler.begin_active();
for (; cell!=endc; ++cell, ++dg_cell)
if (cell->is_locally_owned())
{
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 i=0; i<dim; i++)
for (unsigned int j=0; j<dim; j++)
{
for (unsigned int q=0; q<qf_cell.size(); ++q)
{
local_history_values_at_qpoints[i][j](q)
= lqph[q].get_F()[i][j];
qpoint_to_dof_matrix.vmult (local_history_fe_values[i][j],
local_history_values_at_qpoints[i][j]);
dg_cell->set_dof_values (local_history_fe_values[i][j],
history_field[i][j]);
local_history_values_at_qpoints_stress[i][j](q)
= lqph[q].get_tau()[i][j];
qpoint_to_dof_matrix.vmult (local_history_fe_values_stress[i][j],
local_history_values_at_qpoints_stress[i][j]);
dg_cell->set_dof_values (local_history_fe_values_stress[i][j],
history_field_stress[i][j]);
}
}
}
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(1, DataComponentInterpretation::component_is_scalar);
data_out.add_data_vector(history_dof_handler, history_field[1][1], "F_22",
data_component_interpretation);
data_out.add_data_vector(history_dof_handler, history_field_stress[0][0], "sigma_11",
data_component_interpretation);
data_out.add_data_vector(history_dof_handler, history_field_stress[1][1], "sigma_22",
data_component_interpretation);
data_out.add_data_vector(history_dof_handler, history_field_stress[0][1], "sigma_12",
data_component_interpretation);
//////////////////////////
LA::MPI::Vector soln ;
soln.reinit(locally_owned_dofs,mpi_communicator);
for (unsigned int i = 0; i < soln.size(); ++i)
soln(i) = solution(i);
MappingQEulerian<dim, LA::MPI::Vector > q_mapping(degree, dof_handler, soln);
Vector<float> subdomain (triangulation.n_active_cells());
for (unsigned int i=0; i<subdomain.size(); ++i)
subdomain(i) = triangulation.locally_owned_subdomain();
data_out.add_data_vector (subdomain, "subdomain");
data_out.build_patches(q_mapping, degree);
const unsigned int cycle = time.get_timestep();
const std::string filename = ("solution-" +
Utilities::int_to_string (cycle, 2) +
"." +
Utilities::int_to_string
(triangulation.locally_owned_subdomain(), 4));
std::ofstream output ((filename + ".vtu").c_str());
data_out.write_vtu (output);
if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
{
std::vector<std::string> filenames;
for (unsigned int i=0;
i<Utilities::MPI::n_mpi_processes(mpi_communicator);
++i)
filenames.push_back ("solution-" +
Utilities::int_to_string (cycle, 2) +
"." +
Utilities::int_to_string (i, 4) +
".vtu");
std::ofstream master_output (("solution-" +
Utilities::int_to_string (cycle, 2) +
".pvtu").c_str());
data_out.write_pvtu_record (master_output, filenames);
}
}
}
/////////////////////////////////////////////////////////
int main (int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace PhaseField;
//deallog.depth_console(0);
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
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;
}
parameters.prm
Description: Binary data
