Different boundary constraint behavior in 2D vs 3D (MPI, deal.II 9.4.2)
ME20D503 NEWTON
the solutions for 2d and 3d are attached.
Thank you.
The code works correctly for:
Solid<3>but shows unexpected behavior for:
Solid<2>Expected behavior (both 2D and 3D)Left boundary (boundary_id = 0): fixed in x-direction only
Right boundary (boundary_id = 1): prescribed displacement in x-direction
Other boundaries: free
Entire boundary appears constrained in both x and y
Only the elements directly adjacent to the right boundary show deformation
Interior behaves almost as if fully fixed
The same setup in 3D produces the correct deformation pattern.
please find the code snippet below:
template <int dim>
void Solid<dim>::make_grid()
{
std::vector< unsigned int > repetitions(dim, parameters.elements_per_edge);
repetitions[1] = parameters.elements_per_width;
const Point<dim> bottom_left = (dim == 3 ? Point<dim>(0.0, 0.0, 0.0) : Point<dim>(0.0, 0.0));
const Point<dim> top_right = (dim == 3 ? Point<dim>(parameters.width, parameters.length, parameters.height)
:Point<dim>(parameters.width, parameters.length));
GridGenerator::subdivided_hyper_rectangle(triangulation,
repetitions,
bottom_left,
top_right,
true);
typename Triangulation<dim>::face_iterator
face = triangulation.begin_face(),
endface = triangulation.end_face();
for (; face!=endface; ++face)
if (face->at_boundary())
// if (face->boundary_id() == 0)
{
const Point<dim> face_center (face->center());
if (std::abs(face_center[0]-0.0)<1e-10)
{
std::cout<<"check 0"<<std::endl;
face->set_boundary_id(0);
}
if (std::abs(face_center[0]-parameters.width)<1e-10)
{
std::cout<<"check 1"<<std::endl;
face->set_boundary_id(1);
}
if (std::abs(face_center[1]-0.0)<1e-10)
{
std::cout<<"check 2"<<std::endl;
face->set_boundary_id(2);
}
if (std::abs(face_center[1]-parameters.length)<1e-10)
{
std::cout<<"check 3"<<std::endl;
face->set_boundary_id(3);
}
}
GridTools::scale(parameters.scale, triangulation);
vol_reference = GridTools::volume(triangulation);
vol_current = vol_reference;
pcout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
}
template <int dim>
void Solid<dim>::make_constraints(const int &it_nr)
{
// std::cout << "dim = " << dim << std::endl;
// std::cout << "n_components = " << fe.n_components() << std::endl;
pcout << " CST " << std::flush;
if (it_nr > 1)
return;
constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler_ref, constraints);
const bool apply_dirichlet_bc = (it_nr == 0);
const FEValuesExtractors::Scalar x_displacement(0);
const FEValuesExtractors::Scalar y_displacement(1);
if (time.current() == 0.0)
{
const int boundary_id = 0;
if (apply_dirichlet_bc == true)
{
std::cout<<"bId 3"<<std::endl;
const double length_nm = parameters.length*1.0e-9;
VectorTools::interpolate_boundary_values (dof_handler_ref,
boundary_id,
IncrementalBoundaryValues<dim>(length_nm,
0.0, 0., 0., 0., 0.0, 0., 0., 0., 0.),
constraints,
fe.component_mask(x_displacement));
}
else
{
std::cout<<"bId 4"<<std::endl;
VectorTools::interpolate_boundary_values (dof_handler_ref,
boundary_id,
Functions::ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
}
}
else {
const int boundary_id = 0;
if (apply_dirichlet_bc == true)
{
std::cout<<"bId 3"<<std::endl;
const double length_nm = parameters.length*1.0e-9;
VectorTools::interpolate_boundary_values (dof_handler_ref,
boundary_id,
Functions::ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
}
else
{
std::cout<<"bId 4"<<std::endl;
VectorTools::interpolate_boundary_values (dof_handler_ref,
boundary_id,
Functions::ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
}
}
if (time.current() == 0.0)
{
const int boundary_id = 1;
if (apply_dirichlet_bc == true)
{
std::cout<<"bId 3"<<std::endl;
const double length_nm = parameters.length*1.0e-9;
VectorTools::interpolate_boundary_values (dof_handler_ref,
boundary_id,
IncrementalBoundaryValues<dim>(length_nm,
0.06, 0., 0., 0., 0.0, 0., 0., 0., 0.),
constraints,
fe.component_mask(x_displacement));
}
else
{
std::cout<<"bId 4"<<std::endl;
VectorTools::interpolate_boundary_values (dof_handler_ref,
boundary_id,
Functions::ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
}
}
else {
const int boundary_id = 1;
if (apply_dirichlet_bc == true)
{
std::cout<<"bId 3"<<std::endl;
const double length_nm = parameters.length*1.0e-9;
VectorTools::interpolate_boundary_values (dof_handler_ref,
boundary_id,
Functions::ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
}
else
{
std::cout<<"bId 4"<<std::endl;
VectorTools::interpolate_boundary_values (dof_handler_ref,
boundary_id,
Functions::ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
}
}
constraints.close();
}
template <int dim>
class IncrementalBoundaryValues : public Function<dim>
{
public:
IncrementalBoundaryValues (const double length,
const double applied_strainxx, const double applied_strainxy, const double applied_strainxz,
const double applied_strainyx, const double applied_strainyy, const double applied_strainyz,
const double applied_strainzx, const double applied_strainzy, const double applied_strainzz);
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 double length;
const double applied_strainxx;
const double applied_strainxy;
const double applied_strainxz;
const double applied_strainyx;
const double applied_strainyy;
const double applied_strainyz;
const double applied_strainzx;
const double applied_strainzy;
const double applied_strainzz;
};
template <int dim>
IncrementalBoundaryValues<dim>::
IncrementalBoundaryValues (const double length,
const double applied_strainxx, const double applied_strainxy, const double applied_strainxz,
const double applied_strainyx, const double applied_strainyy, const double applied_strainyz,
const double applied_strainzx, const double applied_strainzy, const double applied_strainzz)
:
Function<dim> (dim),
length (length),
applied_strainxx (applied_strainxx),
applied_strainxy (applied_strainxy),
applied_strainxz (applied_strainxz),
applied_strainyx (applied_strainyx),
applied_strainyy (applied_strainyy),
applied_strainyz (applied_strainyz),
applied_strainzx (applied_strainzx),
applied_strainzy (applied_strainzy),
applied_strainzz (applied_strainzz)
{}
template <int dim>
void
IncrementalBoundaryValues<dim>::
vector_value (const Point<dim> &p,
Vector<double> &values) const
{
Assert (values.size() == dim,
ExcDimensionMismatch (values.size(), dim));
values = 0;
values(0) = applied_strainxx*p[0]+applied_strainxy*p[1];
values(1) = applied_strainyx*p[0]+applied_strainyy*p[1];
// values(0) = applied_strainxx*p[0]+applied_strainxy*p[1]+applied_strainxz*p[2];
// values(1) = applied_strainyx*p[0]+applied_strainyy*p[1]+applied_strainyz*p[2];
// values(2) = applied_strainzx*p[0]+applied_strainzy*p[1]+applied_strainzz*p[2];
}
template <int dim>
void
IncrementalBoundaryValues<dim>::
vector_value_list (const std::vector<Point<dim> > &points,
std::vector<Vector<double> > &value_list) const
{
const unsigned int n_points = points.size();
Assert (value_list.size() == n_points,
ExcDimensionMismatch (value_list.size(), n_points));
for (unsigned int p=0; p<n_points; ++p)
IncrementalBoundaryValues<dim>::vector_value (points[p],
value_list[p]);
}
int main (int argc, char *argv[])
{
using namespace dealii;
using namespace PhaseField;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
try
{
deallog.depth_console(0);
Solid<2> solid_3d("parameters.prm");
solid_3d.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl << exc.what()
<< std::endl << "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl << "Aborting!"
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
