On Mon, May 10, 2010 at 10:44 PM, James Avery <[email protected]> wrote:
> On Thu, May 6, 2010 at 7:22 PM, <[email protected]> wrote:
>> Hi James,
>>
>> actually you'd have to exchange the cases Wolfgang gave below, see the
>> documentation of GeometryInfo<dim>.
>>
>> Furthermore, you might also consider the following alternative approach:
>> a) Compute the volume/area of the cell by numerical quadrature (using
>> FEValues).
>> b) Compute the area/length of the faces of the cell by numerical quadrature
>> (using FEFaceValues).
>> c) Obtain an estimate of the length in one direction by dividing the cell
>> volume by the mean value of the area of the two faces that are orthogonal to
>> that direction (faces 0 and 1 for direction x(or 0)).
>
> Thanks a lot to both of you! I'll have a go at both methods today and
> see what's easiest -- I'll post again with the result.
I've now implemented Wolfgang Bangert's method (with Tobias Leicht's
correction of using the line enumeration described in
GeometryInfo<dim>). I don't know whether I can attach files on the
mailing list, so I'll just copy in the changes instead:
In deal.II/includes/grid/tria_accessor.h, the following has been
inserted in the TriaAccessor class (between the diameter() and
center() declarations):
***************************** snip *****************************
/**
* Length of an object in the direction of
* the given axis, specified in the local
* coordinate system. See the
documentation
* of GeometryInfo<dim> for the meaning
and
* enumeration of the local axes.
*/
double extent_in_direction(const unsigned int axis) const;
***************************** snip *****************************
and in deal.II/sources/grid/tria_accessor.cc:
***************************** snip *****************************
#if deal_II_dimension == 1
template <> double TriaAccessor<1,1,1>::extent_in_direction(const
unsigned int axis) const {
if(axis>0)
Assert (false, ExcIndexRange (direction, 0, 1));
return this->diameter();
}
#elsif deal_II_dimension == 2
template <> double TriaAccessor<2,2,2>::extent_in_direction(const
unsigned int axis) const {
using std::max;
int lines[2][2] = {{2,3}, /// Lines along x-axis -- see
GeometryInfo<dim>, N4), 2D
{0,1}};/// Lines along y-axis
if(axis>1)
Assert (false, ExcIndexRange (direction, 0, 2));
return max(this->line(lines[axis][0])->diameter(),
this->line(lines[axis][1])->diameter());
}
#elsif deal_II_dimension == 3
template <> double TriaAccessor<3,3,3>::extent_in_direction(const
unsigned int axis) const {
using std::max;
int lines[3][4] = {{2,3,6,7}, /// Lines along x-axis -- see
GeometryInfo<dim>, N4), 3D
{0,1,4,5}, /// Lines along y-axis
{8,9,10,11}}; /// Lines along z-axis
double lengths[4];
if(axis>2)
Assert (false, ExcIndexRange (direction, 0, 3));
for(unsigned int i=0;i<4;i++)
lengths[i] = this->line(lines[axis][i])->diameter();
return max(max(lengths[0],lengths[1]), max(lengths[2],lengths[3]));
}
#else
// Not implemented for higher dimensions
#endif
***************************** snip *****************************
I'm not exactly sure how things works with structdim, dim and
spacedim, so I hope I did that part correctly. Anyway, it works quite
well, it seems! I used the following scheme to regularize the aspect
ratios:
***************************** snip *****************************
fespace_member(void) refine_to_aspect_ratio(double max_aspect_ratio)
{
typename DoFHandler<dim>::active_cell_iterator
cell, endc = dof_handler.end();
bool done = false;
while(!done){
done = true;
unsigned int num_refined = 0;
double max_found_ratio = 1.0, Lm[dim];
int im,jm;
for(cell = dof_handler.begin_active();cell!=endc;cell++){
double L[dim];
for(int i=0; i<dim;i++)
L[i] = cell->extent_in_direction(i);
for(int i=0; i<dim;i++)
for(int j=0; j<dim;j++)
if(L[i] > max_aspect_ratio*L[j]){
if(L[i]/L[j] > max_found_ratio){
max_found_ratio = L[i]/L[j];
im = i; jm = j;
for(int k=0;k<dim;k++) Lm[k] = L[k];
}
cell->set_refine_flag(RefinementCase<dim>::cut_axis(i));
done = false;
num_refined++;
}
}
if(dim==3)
fprintf(stderr,"Refining %d cells. Maximal aspect ratio found is %g.
(%g:%g:%g, axis %d>%d)\n",num_refined,max_found_ratio,
Lm[0],Lm[1],Lm[2],im,jm);
triangulation.execute_coarsening_and_refinement();
}
}
***************************** snip *****************************
An example of the result with max_aspect_ratio=3 can be seen at
http://diku.dk/~avery/aspect-out.tar.gz, where the important step is
between initial-mesh.msh and regularized-mesh.msh.
Thanks a lot for the help, and I hope the additions to
tria_accessor.{h,cpp} are useful!
--
Med venlig hilsen,
James Avery <[email protected]>
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
