Thank you very much for your help. After some trial and error, I was able
to do it. I have attached the code if that can be of any help to others
searching for similar solutions. I was modifying step-1.cc and did not
modify the file name. Thank you for your suggestion of using the error
estimation for corsening, I will try that too.
Kaushik
On Mon, Nov 30, 2020 at 2:52 PM Wolfgang Bangerth <[email protected]>
wrote:
> On 11/30/20 12:36 PM, Kaushik Das wrote:
> > saveRefinement<3> SaveRefinement(triangulation);
>
> Unrelated to your question, but this convention is very hard to read for
> most
> programmers: You name a class with lower-case first letter and a variable
> with
> upper-case first letter. Both use CamelCase. It's quite difficult to see
> from
> the name what's a variable and what's a class.
>
>
> > triangulation.signals.pre_refinement.connect([&SaveRefinement] {
> > SaveRefinement(); })
> >
> > /// refine the mesh
> > ....
> > triangulation.load_coarsen_flags(SaveRefinement.cellRefined);
> > triangulation.execute_coarsening_and_refinement();
> >
> > At this point, I am hitting Assert(v.size() == n_active_cells(),
> > ExcGridReadError()) in load_coarsen_flags.
> >
> > Which also makes sense to me now, because those flags saved from
> > save_refine_flags are now parent cells, and I have to mark their child
> cells
> > for coarsening.
>
> That's correct.
>
>
> > I will greatly appreciate any hint on how that can be done.
>
> You need to store a map from the old cell to the flags. You could use
> std::map<Triangulation::cell_iterator,bool>
> for example.
>
>
> > I looked into the Step-26 example. There is another reason I wanted to
> refine
> > the mesh based on geometric criteria. In addition to the laser spot, I
> would
> > like to refine the mesh near the layer that is being printed at the
> given time
> > step. And coarsen the layers underneath it. I am trying to simulation a
> layer
> > upon layer printing process following a paper by Claire Bruna-Rosso et
> al.,
> > who also used deal.ii. They used the FE_Nothing element to exclude the
> layers
> > that are yet to be printed from the solution.
>
> You can refine by your geometric argument, but it might be easiest to
> coarsen
> based on the general indicators.
>
> Best
> W.
>
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> www: http://www.math.colostate.edu/~bangerth/
>
> --
> 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 a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/6hhi5Xf0_mQ/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> [email protected].
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/696423e9-f1ac-6394-04a4-676fb6668040%40colostate.edu
> .
>
--
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/CAC-fs6sOCQ%3DOOXAhbcp4prvPQvWzBKUDRVqf8ZRHW0oUPwGycg%40mail.gmail.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 1999 - 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.
*
* ---------------------------------------------------------------------
*/
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <iostream>
#include <fstream>
#include <cmath>
using namespace dealii;
template <int dim>
void vtk_out(const Triangulation<dim> &triangulation, int level = -1)
{
std::string fileName =
"grid-" + (level >= 0 ? std::to_string(level) : "") + ".vtk";
std::ofstream out("grid-" + std::to_string(level) + ".vtk");
GridOut grid_out;
grid_out.write_vtk(triangulation, out);
std::cout << "Grid written to " + fileName << std::endl;
}
template <int dim, int spacedim = dim>
class SaveRefinedCells
{
private:
std::vector<std::vector<TriaIterator<CellAccessor<dim, spacedim>>>>
cells_refined_at_lavel;
unsigned int refinement_level;
Triangulation<dim, spacedim> &triangulation;
void set_coarsen_flag_to_child(TriaIterator<CellAccessor<dim, spacedim>> cell)
{
if (!cell->is_active())
{
for (unsigned int c = 0; c < cell->n_children(); c++)
{
set_coarsen_flag_to_child(cell->child(c));
}
}
else
{
cell->set_coarsen_flag();
}
}
public:
SaveRefinedCells(Triangulation<dim, spacedim> &tri)
: triangulation(tri)
, refinement_level(0)
{}
void operator()()
{
for (auto cell : triangulation.cell_iterators())
{
if (cell->refine_flag_set())
{
cells_refined_at_lavel[refinement_level].push_back(cell);
}
}
}
void setRefinementLevel(unsigned int L)
{
refinement_level = L;
if (refinement_level + 1 > cells_refined_at_lavel.size())
{
cells_refined_at_lavel.resize(refinement_level + 1);
}
}
auto &getCellsRefined()
{
return cells_refined_at_lavel;
}
void clear()
{
cells_refined_at_lavel.clear();
}
void revert_refinement()
{
auto &celllists = getCellsRefined();
for (auto i = celllists.rbegin(); i != celllists.rend(); ++i)
{
auto cells = *i;
for (auto cell : cells)
{
set_coarsen_flag_to_child(cell);
}
triangulation.execute_coarsening_and_refinement();
}
}
};
template <int dim, int spacedim = dim>
double refine_cylinder(Triangulation<dim, spacedim> &triangulation,
const Point<dim> & cylCenter,
const double height,
const double radius,
const double initElementSize)
{
double endElementSize = initElementSize;
constexpr unsigned int maxRefinementSteps(4);
for (unsigned int step = 0; step < maxRefinementSteps; ++step)
{
double disCriteria = 1.001 * (endElementSize / 2) * std::sqrt(2);
for (auto cell : triangulation.active_cell_iterators())
{
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
++v)
{
const auto &vertex = cell->vertex(v);
Point<3> center(cylCenter(0), cylCenter(1), vertex(2));
const double distance_from_center = center.distance(vertex);
if (std::fabs(distance_from_center - radius) <= disCriteria)
{
cell->set_refine_flag();
break;
}
}
}
triangulation.execute_coarsening_and_refinement();
endElementSize /= 2;
}
return endElementSize;
}
template <int dim, int spacedim = dim>
void refine_near_layer(Triangulation<dim, spacedim> &triangulation,
const double layerZ,
const double initElementSize,
SaveRefinedCells<dim> & save_refined_cells)
{
double elementSize = initElementSize;
for (unsigned int step = 0; step < 2; ++step)
{
save_refined_cells.setRefinementLevel(step);
double dist = 1.001 * (elementSize / 2) * std::sqrt(2);
for (auto cell : triangulation.active_cell_iterators())
{
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
++v)
{
const auto & vertex = cell->vertex(v);
const double distance_from_layer = std::abs(vertex(2) - layerZ);
if (distance_from_layer <= dist)
{
cell->set_refine_flag();
break;
}
}
}
triangulation.execute_coarsening_and_refinement();
elementSize /= 2;
}
}
int main()
{
Triangulation<3> triangulation;
double height{10};
double initElementSize{2};
double radius{1};
double buffer{1};
double length = buffer + 2 * radius;
Point<3> cylCenter{length / 2, length / 2, 0};
Point<3> left{0, 0, 0};
Point<3> right{length, length, height};
std::vector<unsigned int> repetitions{
(unsigned int)(length / initElementSize),
(unsigned int)(length / initElementSize),
(unsigned int)(height / initElementSize)};
GridGenerator::subdivided_hyper_rectangle(
triangulation, repetitions, left, right, true);
//
// Initial refinement
//
auto element_size_after_init_refinement =
refine_cylinder(triangulation, cylCenter, height, radius, initElementSize);
//
// Attach a signal that will save the cells for refinement in subsequest
// refinments
//
SaveRefinedCells<3> save_refined_cells(triangulation);
triangulation.signals.pre_refinement.connect(
[&save_refined_cells] { save_refined_cells(); });
vtk_out(triangulation, 0);
double layerZ = 2;
for (int i = 0; i < 5; i++)
{
double elemSize = element_size_after_init_refinement;
//
// Refine again at the start of the loop
//
refine_near_layer(triangulation, layerZ, elemSize, save_refined_cells);
vtk_out(triangulation, i + 1);
//
// revert the refinement done at the start of the loop
//
save_refined_cells.revert_refinement();
save_refined_cells.clear();
layerZ += 2;
}
vtk_out(triangulation, 6);
}
