1 #include "metric_to_mesh_generator.h" 2 #include <deal.II/dofs/dof_renumbering.h> 3 #include <deal.II/distributed/shared_tria.h> 4 #include <deal.II/dofs/dof_tools.h> 5 #include <deal.II/lac/la_parallel_vector.h> 8 #include "grid_refinement/gmsh_out.h" 13 template<
int dim,
int nstate,
typename real,
typename MeshType>
15 std::shared_ptr<dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType>> _volume_nodes_mapping,
16 std::shared_ptr<MeshType> _triangulation)
17 : volume_nodes_mapping(_volume_nodes_mapping)
18 , triangulation(_triangulation)
19 , dof_handler_vertices(*triangulation)
21 , fe_system(dealii::FESystem<dim>(fe_q,1))
22 , mpi_communicator(MPI_COMM_WORLD)
23 , pcout(
std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0)
24 , filename(
"mesh_to_be_generated")
25 , filename_pos(filename +
".pos")
26 , filename_geo(filename +
".geo")
27 , filename_msh(filename +
".msh")
34 pcout<<
"This class is currently hardcoded for dim = 2. To be updated later when required."<<std::endl<<std::flush;
41 template<
int dim,
int nstate,
typename real,
typename MeshType>
55 for(
unsigned int i=0; i<dim*dim; ++i)
64 template<
int dim,
int nstate,
typename real,
typename MeshType>
66 const std::vector<dealii::Tensor<2, dim, real>> &cellwise_optimal_metric)
69 AssertDimension(cellwise_optimal_metric.size(),
triangulation->n_active_cells());
71 dealii::Quadrature<dim> quadrature =
fe_system.get_unit_support_points();
73 const unsigned int n_dofs_cell =
fe_system.dofs_per_cell;
74 std::vector<dealii::types::global_dof_index> dof_indices(n_dofs_cell);
75 AssertDimension(n_dofs_cell, dealii::GeometryInfo<dim>::vertices_per_cell)
80 if(! cell->is_locally_owned()) {
continue;}
82 fe_values_vertices.reinit(cell);
83 cell->get_dof_indices(dof_indices);
85 const std::vector<dealii::Point<dim>> &cell_vertices = fe_values_vertices.get_quadrature_points();
87 for(
unsigned int idof = 0; idof < cell_vertices.size(); ++idof)
89 const unsigned int iquad = idof;
93 AssertDimension(cell_vertices.size(), dealii::GeometryInfo<dim>::vertices_per_cell);
96 for(
unsigned int idof = 0; idof < cell_vertices.size(); ++idof)
98 const dealii::Point<dim> &ref_vertex = dealii::GeometryInfo<dim>::unit_cell_vertex(idof);
99 const dealii::Point<dim> expected_vertex =
volume_nodes_mapping->transform_unit_to_real_cell(cell, ref_vertex);
100 const real error_in_vertex = expected_vertex.distance(
all_vertices[dof_indices[idof]]);
101 if(error_in_vertex > 1.0e-10)
103 std::cout<<
"Expected vertex = "<<expected_vertex<<
" "<<
"Actual vertex = "<<
all_vertices[dof_indices[idof]]<<std::endl;
104 std::cout<<
"Error in ordering of vertices. Aborting.."<<std::flush;
113 metric_count_at_vertices = 0;
117 if(! cell->is_locally_owned()) {
continue;}
119 cell->get_dof_indices(dof_indices);
120 const unsigned int cell_index = cell->active_cell_index();
122 for(
unsigned int idof = 0; idof<n_dofs_cell; ++idof)
124 const unsigned int idof_global = dof_indices[idof];
125 metric_count_at_vertices(idof_global) += 1;
126 for(
unsigned int i=0; i<dim; ++i)
128 for(
unsigned int j=0; j<dim; ++j)
137 metric_count_at_vertices.compress(dealii::VectorOperation::add);
138 metric_count_at_vertices.update_ghost_values();;
139 for(
unsigned int i=0; i<dim*dim; ++i)
148 for(
unsigned int i=0; i<dim*dim; ++i)
155 for(
unsigned int i=0; i<dim*dim; ++i)
161 template<
int dim,
int nstate,
typename real,
typename MeshType>
164 AssertDimension(
fe_system.dofs_per_cell, dealii::GeometryInfo<dim>::vertices_per_cell);
169 const std::string quotes =
"\"";
170 const std::string filename_this_processor =
filename + std::to_string(
mpi_rank) +
".pos";
171 std::ofstream outfile(filename_this_processor);
184 std::vector<dealii::types::global_dof_index> dof_indices(
fe_system.dofs_per_cell);
188 if(! cell->is_locally_owned()) {
continue;}
189 cell->get_dof_indices(dof_indices);
191 for(
const auto &tri :
192 std::array<std::array<int,3>,2>{{
198 for(
unsigned int i = 0; i < tri.size(); ++i)
200 if(i != 0) {outfile <<
",";}
202 const unsigned int idof = tri[i];
203 const unsigned int idof_global = dof_indices[idof];
205 const dealii::Point<dim> &pos =
all_vertices[idof_global];
206 if(dim >= 1){outfile << pos[0] <<
",";}
207 else {outfile << 0 <<
",";}
209 if(dim >= 2){outfile << pos[1] <<
",";}
210 else {outfile << 0 <<
",";}
212 if(dim >= 3){outfile << pos[2];}
219 for(
unsigned int i = 0; i < tri.size(); ++i)
221 const unsigned int idof = tri[i];
222 const unsigned int idof_global = dof_indices[idof];
224 dealii::Tensor<2,dim,real> vertex_metric;
225 for(
unsigned int idim =0; idim < dim; ++idim)
227 for(
unsigned int jdim =0; jdim < dim; ++jdim)
233 for(
unsigned int idim =0; idim < 3; ++idim)
235 for(
unsigned int jdim =0; jdim < 3; ++jdim)
237 if(flag) {outfile<<
",";}
240 if( (idim < dim) && (jdim < dim) )
242 outfile<<vertex_metric[idim][jdim];
247 if(idim == jdim) {outfile<<1.0;}
267 outfile0<<
"// Background mesh containing optimal metric field at nodes."<<
'\n';
268 outfile0<<
"View "<<quotes<<
"nodalMetric"<<quotes<<
" {"<<
'\n';
270 for(
int iproc = 0; iproc <
n_mpi; ++iproc)
272 const std::string input_filename =
filename + std::to_string(iproc) +
".pos";
273 std::ifstream inputfile (input_filename);
274 outfile0<<inputfile.rdbuf();
277 outfile0<<
"};"<<
'\n';
285 template<
int dim,
int nstate,
typename real,
typename MeshType>
291 outfile<<
" // GEO file "<<
'\n'<<
'\n';
292 const std::string quotes =
"\"";
294 outfile<<
"SetFactory("<<quotes<<
"OpenCASCADE"<<quotes<<
");"<<
'\n'<<
'\n';
300 outfile<<
"Merge "<<quotes<<
filename_pos<<quotes<<
";"<<
'\n';
303 outfile<<
"Background Mesh View[0];"<<
'\n';
306 outfile<<
"Mesh.SmoothRatio = 0;"<<
'\n';
307 outfile<<
"Mesh.AnisoMax = 1e30;"<<
'\n';
308 outfile<<
"Mesh.Algorithm = 7;"<<
'\n';
311 outfile<<
"Mesh.RecombinationAlgorithm = 2;"<<
'\n';
312 outfile<<
"Mesh.RecombineAll = 1;"<<
'\n';
317 template<
int dim,
int nstate,
typename real,
typename MeshType>
319 const std::vector<dealii::Tensor<2, dim, real>> &cellwise_optimal_metric)
330 int val = std::system(command_str.c_str());
332 pcout<<
"Generated new mesh."<<std::endl;
337 template<
int dim,
int nstate,
typename real,
typename MeshType>
const std::string filename_msh
.geo file
dealii::DoFHandler< dim > dof_handler_vertices
DoFHandler to get global index of vertices.
MetricToMeshGenerator(std::shared_ptr< dealii::MappingFEField< dim, dim, VectorType, DoFHandlerType >> _volume_nodes_mapping, std::shared_ptr< MeshType > _triangulation)
Constructor.
dealii::IndexSet ghost_vertex_dofs
Dofs not owned by current processor (but it can still access them).
std::shared_ptr< MeshType > triangulation
Stores triangulation currently in use.
const dealii::FESystem< dim > fe_system
FESystem for vertices, created with nstate = 1 to relate an entire vertex of size dim by a single dof...
Files for the baseline physics.
dealii::IndexSet locally_owned_vertex_dofs
Dofs owned by current processor.
int mpi_rank
Processor# of current processor.
static void write_geo_hyper_cube(const double left, const double right, std::ostream &out, const bool colorize=true)
Writes the part of the .geo file associated wit the hyperdube geometry.
void generate_mesh_from_cellwise_metric(const std::vector< dealii::Tensor< 2, dim, real >> &cellwise_optimal_metric)
Generates mesh from the input cellwise metric.
std::shared_ptr< dealii::MappingFEField< dim, dim, VectorType, DoFHandlerType > > volume_nodes_mapping
Mapping field to update physical quadrature points, jacobians etc with the movement of volume nodes...
void write_pos_file() const
Writes .pos file.
void reinit()
Reinitialize dof handler vertices after updating triangulation.
dealii::IndexSet locally_relevant_vertex_dofs
Union of locally owned degrees of freedom and ghost degrees of freedom.
std::vector< dealii::Point< dim > > all_vertices
Stores all vertices.
void interpolate_metric_to_vertices(const std::vector< dealii::Tensor< 2, dim, real >> &cellwise_optimal_metric)
Interpolates cellwise metric field to vertices.
void write_geo_file() const
Writes .geo file.
std::string get_generated_mesh_filename() const
Returns name of the .msh file.
int n_mpi
Total no. of processors.
std::array< VectorType, dim *dim > optimal_metric_at_vertices
Stores optimal metric at vertices. Accessed by optimal_metric_at_vertices[position_in_metric_tensor][...
const std::string filename_pos
.pos file.
dealii::ConditionalOStream pcout
std::cout only on processor #0.
MPI_Comm mpi_communicator
Alias for MPI_COMM_WORLD.
Class to convert metric field to mesh using BAMG.
const std::string filename_geo
.geo file
const std::string filename
Name of the file.