Hello everyone, I have been reading the messages in this user group for a while, but this is my first post! I am hoping that you can help with my problem.
I am coming from a background of solid mechanics FEA and I've been using
CalculiX as an openSource program for some time. My group has decided to
switch to deal.ii to implement some advanced functionality that the other
programs don't have. I'm an engineer with some programming experience,
though fairly new to C++.
We are working with some complex geometries that are not practical to mesh
with only hex elements. The meshes are generated externally, so I am
trying to import the meshes into deal.ii. I have tried to do this in two
general ways and run into problems with each, so I will split my post into
2 parts.
Part 1: Abaqus .inp format
I had success importing hex meshes from the abaqus <.inp> format, but it
seems that the UCD reader does not accept tet meshes. I am basing this on
the following bit of code in grid_in::read_ucd around Line 948
if (((cell_type == "line") && (dim == 1)) ||
((cell_type == "quad") && (dim == 2)) ||
((cell_type == "hex") && (dim == 3)))
// found a cell
{
The absence of
((cell_type == "tet") && (dim == 3))
suggests that simplex elements are not supported when importing from .inp.
Am I correct in this conclusion. How difficult would it be for me to add
this functionality in?
Part 2: GMSH format
While looking at grid_in.cc, I saw that the MSH reader supported tet
elements. I found the FAQ page supplied these examples of simplex meshes:
https://www.dealii.org/developer/doxygen/deal.II/group__simplex.html
I first ran the 2D simplex mesh example, and then I modified it to work for
3D without much problem. Then, I switched to the mixed mesh, and had
success running the 2D example. However, I started having problems when I
switched to 3D. The mesh seems to import, but then I have issues operating
with it. I found several places within fe_simplex_p.cc (L714, L759) where
the dimension is asserted to be 2:
AssertDimension(dim, 2);
This seems like a problem for a 3D analysis, so I am wondering if I need to
do something different so that this function doesn't get called.
The problem seems to arrive when I call dof_handler.distribute_dofs(fe);
I've attached my modified version of the example program along with my
simple mixed-element 3D mesh.
Is there something that I am missing regarding the dof_handler when trying
to use a 3D mixed-element mesh?
Thanks very much,
Alex
--
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/a905250a-c35f-4a75-952a-a2237451e94fn%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 1999 - 2021 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.
*
* ---------------------------------------------------------------------
*/
#include <deal.II/base/function.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <fstream>
#include <iostream>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/hp/q_collection.h>
using namespace dealii;
class Step3
{
public:
Step3();
void run();
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
Triangulation<3> triangulation;
const hp::MappingCollection<3> mapping;
const hp::FECollection<3> fe;
const hp::QCollection<3> quadrature_formula;
DoFHandler<3> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
Step3::Step3()
: mapping(MappingFE<3>(FE_SimplexP<3>(1)), MappingFE<3>(FE_Q<3>(1)))
, fe(FE_SimplexP<3>(2), FE_Q<3>(2))
, quadrature_formula(QGaussSimplex<3>(3), QGauss<3>(3))
, dof_handler(triangulation)
{}
void Step3::make_grid()
{
GridIn<3>(triangulation).read("input/mixed.msh");
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
void Step3::setup_system()
{
for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->reference_cell() == ReferenceCells::Tetrahedron)
cell->set_active_fe_index(0);
else if (cell->reference_cell() == ReferenceCells::Hexahedron)
cell->set_active_fe_index(1);
else
Assert(false, ExcNotImplemented());
}
dof_handler.distribute_dofs(fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
void Step3::assemble_system()
{
hp::FEValues<3> hp_fe_values(mapping,
fe,
quadrature_formula,
update_values | update_gradients |
update_JxW_values);
FullMatrix<double> cell_matrix;
Vector<double> cell_rhs;
std::vector<types::global_dof_index> local_dof_indices;
for (const auto &cell : dof_handler.active_cell_iterators())
{
hp_fe_values.reinit(cell);
const auto &fe_values = hp_fe_values.get_present_fe_values();
const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
cell_rhs.reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
cell_matrix = 0;
cell_rhs = 0;
for (const unsigned int q_index : fe_values.quadrature_point_indices())
{
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
cell_matrix(i, j) +=
(fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
for (const unsigned int i : fe_values.dof_indices())
cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
1. * // f(x_q)
fe_values.JxW(q_index)); // dx
}
cell->get_dof_indices(local_dof_indices);
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
for (const unsigned int i : fe_values.dof_indices())
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(
mapping, dof_handler, 0, Functions::ZeroFunction<3>(), boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
void Step3::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}
void Step3::output_results() const
{
DataOut<3> data_out;
DataOutBase::VtkFlags flags;
flags.write_higher_order_cells = true;
data_out.set_flags(flags);
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.build_patches(mapping, 2);
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
void Step3::run()
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
int main()
{
deallog.depth_console(2);
Step3 laplace_problem;
laplace_problem.run();
return 0;
}
mixed.msh
Description: Mesh model
