Hi Dealii Users; I have a similar problem, I am running a 3D problem for a beam with the following dimensions x=10 mm and y= 1 mm and z =1 mm I am trying to adapt step-8 to my case the I am facing the following challenges :
1- importing msh file to the dealii, I copied a print_mesh_info function to
see the vtu file but the following error pop up.
An error occurred in line <2018> of file
</home/mohammed/FEM/dealii/dealii/source/grid/grid_in.cc> in function
void dealii::GridIn<dim, spacedim>::read_msh(std::istream&) [with int
dim = 3; int spacedim = 3; std::istream = std::basic_istream<char>]
The violated condition was:
cells.size() > 0
2- the purpose of this work is to learn how to implement the boundary
conditions. the .msh file has been meshed as hexahedral where it has 100
cells along x axis and 10 cells along y and 10 cells along z axis. I would
like to apply 100 N force on the right side of the beam pointing to ( + x
axis and the second case + y axis " I mean the direction of the force") and
a fixed wall on the left side knowing that I named the cells for the right
side of the beam (in .msh file) force while for the left side of the beam
wall.
Thank you in advance !
Kind Regards.
On Tuesday, December 10, 2019 at 4:35:20 PM UTC+1, Hasan Tasneem wrote:
>
> I am a beginner on Deal.ii. I have installed and run some of the
> tutorials. I am now trying to run one specific problem. This problem is 2D
> and is a beam problem. My dimension of the beam is 3*100mm. The beam is
> fixed on one end and there is 1000N of transverse loading on the other end.
> I am following Step-8 which solves an elastic equation. Everything I have
> inputted correctly. However, I am stuck on applying the fixed support. I
> want the fixed support at x equal to 0. I don't know how to proceed
> further.
>
> VectorTools::interpolate_boundary_values
> <https://dealii.org/developer/doxygen/deal.II/namespaceVectorTools.html#a187aeb575be07bc47cb3dea1a47aaf88>
> (dof_handler,
> 0,
> Functions::ZeroFunction<dim>
> <https://dealii.org/developer/doxygen/deal.II/classFunctions_1_1ZeroFunction.html>
> (dim),
> boundary_values);
> MatrixTools::apply_boundary_values
> <https://dealii.org/developer/doxygen/deal.II/namespaceMatrixTools.html#a9ad0eb7a8662628534586716748d62fb>
> (boundary_values,
> system_matrix,
> solution,
> system_rhs);
>
> How should I alter this to get fixed support at x equal to 0.
>
>
>
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/db40fc39-9631-4eed-90e8-0def8efb1441%40googlegroups.com.
##
# CMake script for the step-8 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "beam")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
FIND_PACKAGE(deal.II 9.2.0 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
beam.msh
Description: Mesh model
/* ---------------------------------------------------------------------
*
* Copyright (C) 2000 - 2019 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, University of Heidelberg, 2000
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/tensor.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
// Grid Includes
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <fstream>
#include <iostream>
namespace beam
{
using namespace dealii;
template <int dim>
class ElasticProblem
{
public:
ElasticProblem();
~ElasticProblem();
void run();
private:
void setup_system();
void assemble_system();
void solve();
void refine_grid();
void output_results(const unsigned int cycle) const;
Triangulation<dim> beam;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
AffineConstraints<double> hanging_node_constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
void right_hand_side(const std::vector<Point<dim>> &points,
std::vector<Tensor<1, dim>> & values)
{
Assert(values.size() == points.size(),
ExcDimensionMismatch(values.size(), points.size()));
Assert(dim >= 2, ExcNotImplemented());
Point<dim> point_1, point_2;
point_1(0) = 0.5;
point_2(0) = -0.5;
for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
{
if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) ||
((points[point_n] - point_2).norm_square() < 0.2 * 0.2))
values[point_n][0] = 1.0;
else
values[point_n][0] = 0.0;
if (points[point_n].norm_square() < 0.2 * 0.2)
values[point_n][1] = 1.0;
else
values[point_n][1] = 0.0;
}
}
template <int dim>
ElasticProblem<dim>::ElasticProblem()
: dof_handler(beam)
, fe(FE_Q<dim>(1), dim)
{}
template <int dim>
ElasticProblem<dim>::~ElasticProblem()
{
dof_handler.clear();
}
template <int dim>
void ElasticProblem<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
hanging_node_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler,
hanging_node_constraints);
hanging_node_constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
hanging_node_constraints,
/*keep_constrained_dofs = */ true);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void ElasticProblem<dim>::assemble_system()
{
QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
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> lambda_values(n_q_points);
std::vector<double> mu_values(n_q_points);
Functions::ConstantFunction<dim> lambda(1.), mu(1.);
std::vector<Tensor<1, dim>> rhs_values(n_q_points);
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit(cell);
lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
mu.value_list(fe_values.get_quadrature_points(), mu_values);
right_hand_side(fe_values.get_quadrature_points(), rhs_values);
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;
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
cell_matrix(i, j) +=
( //
(fe_values.shape_grad(i, q_point)[component_i] * //
fe_values.shape_grad(j, q_point)[component_j] * //
lambda_values[q_point]) //
+ //
(fe_values.shape_grad(i, q_point)[component_j] * //
fe_values.shape_grad(j, q_point)[component_i] * //
mu_values[q_point]) //
+ //
((component_i == component_j) ? //
(fe_values.shape_grad(i, q_point) * //
fe_values.shape_grad(j, q_point) * //
mu_values[q_point]) : //
0) //
) * //
fe_values.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;
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
cell_rhs(i) += fe_values.shape_value(i, q_point) *
rhs_values[q_point][component_i] *
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
hanging_node_constraints.condense(system_matrix);
hanging_node_constraints.condense(system_rhs);
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<dim>(dim),
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
template <int dim>
void ElasticProblem<dim>::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<> cg(solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
cg.solve(system_matrix, solution, system_rhs, preconditioner);
hanging_node_constraints.distribute(solution);
}
template <int dim>
void ElasticProblem<dim>::refine_grid()
{
Vector<float> estimated_error_per_cell(beam.n_active_cells());
KellyErrorEstimator<dim>::estimate(
dof_handler,
QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_number(beam,
estimated_error_per_cell,
0.3,
0.03);
beam.execute_coarsening_and_refinement();
}
template <int dim>
void ElasticProblem<dim>::output_results(const unsigned int cycle) const
{
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
std::vector<std::string> solution_names;
switch (dim)
{
case 1:
solution_names.emplace_back("displacement");
break;
case 2:
solution_names.emplace_back("x_displacement");
solution_names.emplace_back("y_displacement");
break;
case 3:
solution_names.emplace_back("x_displacement");
solution_names.emplace_back("y_displacement");
solution_names.emplace_back("z_displacement");
break;
default:
Assert(false, ExcNotImplemented());
}
data_out.add_data_vector(solution, solution_names);
data_out.build_patches();
std::ofstream output("solution-" + std::to_string(cycle) + ".vtk");
data_out.write_vtk(output);
}
template <int dim>
void print_mesh_info(const Triangulation<dim> &triangulation,
const std::string & filename)
{
std::cout << "Mesh info:" << std::endl
<< " dimension: " << dim << std::endl
<< " no. of cells: " << triangulation.n_active_cells() << std::endl;
{
std::map<types::boundary_id, unsigned int> boundary_count;
for (const auto &cell : triangulation.active_cell_iterators())
{
for (const auto &face : cell->face_iterators())
{
if (face->at_boundary())
boundary_count[face->boundary_id()]++;
}
}
std::cout << " boundary indicators: ";
for (const std::pair<const types::boundary_id, unsigned int> &pair :
boundary_count)
{
std::cout << pair.first << "(" << pair.second << " times) ";
}
std::cout << std::endl;
}
std::ofstream out(filename);
GridOut grid_out;
grid_out.write_vtu(triangulation, out);
std::cout << " written to " << filename << std::endl << std::endl;
}
template <int dim>
void ElasticProblem<dim>::run()
{
for (unsigned int cycle = 0; cycle < 8; ++cycle)
{
std::cout << "Cycle " << cycle << ':' << std::endl;
if (cycle == 0)
{
Triangulation<3> beam;
GridIn<3> gridin;
gridin.attach_triangulation(beam);
std::ifstream B("beam.msh");
gridin.read_msh(B);
print_mesh_info(beam,"beam.vtu");
}
else
refine_grid();
std::cout << " Number of active cells: "
<< beam.n_active_cells() << std::endl;
setup_system();
std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
assemble_system();
solve();
output_results(cycle);
}
}
} // namespace Step8
int main()
{
try
{
beam::ElasticProblem<3> elastic_problem_3d;
elastic_problem_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;
}
