Also, here is the new program, with the added if and
prepare_coarsening_and_refinement(). You will find the interesting part
between lines 85 and 104.
On Monday, 18 September 2017 16:29:09 UTC+2, Lucas Campos wrote:
>
> Dear Bruno,
>
> You will find attached the resulting grids. While I originally expected to
> find new divisions only along the x-direction, your previous answer tells
> me it actually cuts in some local coordinate (whatever that might be). In
> this case, what is the best way to refine the cylinder in a single
> direction? I tried adding an if(cell->boundary_id() == 0) before setting
> the flag, but it did not really help. The mesh is still being refined in
> all directions.
>
> I ask this because I plan to do some transformations on this cylinder, to
> turn it into a coil. In this case, the final result is much better if there
> are enough cuts in the x-axis.
>
> Bests
>
> On Monday, 18 September 2017 14:44:40 UTC+2, Bruno Turcksin wrote:
>>
>> Lucas,
>>
>> On Monday, September 18, 2017 at 3:50:20 AM UTC-4, Lucas Campos wrote:
>>>
>>> I am trying to refine a cylinder in the x direction only, but it seems
>>> that when I use
>>>
>>> > cell->set_refine_flag(RefinementCase<3>::cut_x);
>>>
>>> the whole cylinder is refined. Indeed, if I run the above line for every
>>> active cell *once* I would expect the number of cells to double. However,
>>> they quintuple!
>>>
>> It looks good but how does the mesh look like? Here
>> <http://dealii.org/developer/doxygen/deal.II/classCellAccessor.html#afb6cc537720a5b6381c237abe0887de2>,
>>
>> you can see that cut_x is in the local coordinate system not the global
>> one. So I would not be surprise if the mesh doesn't look like you think it
>> should. Also if you do more than one refinement you should also use
>> prepare_coarsening_and_refinement() before you refine your mesh.
>>
>> Best,
>>
>> Bruno
>>
>
--
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].
For more options, visit https://groups.google.com/d/optout.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2013 - 2017 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 at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*
* Author: Timo Heister, Texas A&M University, 2013
*/
#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_tools.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <iostream>
#include <fstream>
#include <map>
using namespace dealii;
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;
// Next loop over all faces of all cells and find how often each
// boundary indicator is used (recall that if you access an element
// of a std::map object that doesn't exist, it is implicitly created
// and default initialized -- to zero, in the current case -- before
// we then increment it):
{
std::map<unsigned int, unsigned int> boundary_count;
typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
{
if (cell->face(face)->at_boundary())
boundary_count[cell->face(face)->boundary_id()]++;
}
}
std::cout << " boundary indicators: ";
for (std::map<unsigned int, unsigned int>::iterator it=boundary_count.begin();
it!=boundary_count.end();
++it)
{
std::cout << it->first << "(" << it->second << " times) ";
}
std::cout << std::endl;
}
// Finally, produce a graphical representation of the mesh to an output
// file:
{
std::ofstream out (filename + ".vtk");
GridOut grid_out;
grid_out.write_vtk (triangulation, out);
}
std::cout << " written to " << filename
<< std::endl
<< std::endl;
}
void grid_cylinder()
{
Triangulation<3> triangulation;
static const CylindricalManifold<3> boundary;
GridGenerator::cylinder (triangulation, 0.05, 2);
triangulation.set_manifold (0, boundary);
triangulation.prepare_coarsening_and_refinement();
print_mesh_info (triangulation, "pre_grid");
for (int i=0; i<1; i++) {
for(auto cell: triangulation.active_cell_iterators()) {
if(cell->boundary_id() == 0)
cell->set_refine_flag(RefinementCase<3>::cut_x);
}
triangulation.execute_coarsening_and_refinement();
}
print_mesh_info (triangulation, "grid");
}
int main ()
{
try
{
grid_cylinder ();
}
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;
}
}
