// not active. In this case, simply
// return the subface_no.
return subface_no;
-
+
const bool first_child_has_children=face->child(0)->has_children();
// if the first child has children
// (FaceRefineCase case_x1y or case_y1x),
}
-
+
// given the number of face's child
// (subface_no) and grandchild
// (subsubface_no), return the number of the
// This can only be an anisotropic refinement case
Assert(face->refinement_case() < RefinementCase<2>::isotropic_refinement,
ExcInternalError());
-
+
const bool first_child_has_children=face->child(0)->has_children();
static const unsigned int e = deal_II_numbers::invalid_unsigned_int;
-
+
// array containing the translation of the
// numbers,
//
{e,1}}, // first subface, second subsubface, first_child_has_children==no and yes
{{1,2}, // second subface, first subsubface, first_child_has_children==no and yes
{2,3}}}; // second subface, second subsubface, first_child_has_children==no and yes
-
+
Assert(translated_subface_no[subface_no][subsubface_no][first_child_has_children]!=e,
ExcInternalError());
= 1/|K| \int_{\hat K} \vec x(xi,eta) |det J| d(xi) d(eta)
# x and y are arrays holding the x- and y-values of the four vertices
- # of this cell in real space.
+ # of this cell in real space.
x := array(0..3);
y := array(0..3);
tphi[0] := (1-xi)*(1-eta):
detJ should be the same, as a matrix and the matrix created from it
by exchanging two consecutive lines and two neighboring columns have
the same determinant.
-
+
# x, y and z are arrays holding the x-, y- and z-values of the four vertices
- # of this cell in real space.
+ # of this cell in real space.
x := array(0..7):
y := array(0..7):
z := array(0..7):
with (linalg):
J := matrix(3,3, [[diff(x_real, xi), diff(x_real, eta), diff(x_real, zeta)],
[diff(y_real, xi), diff(y_real, eta), diff(y_real, zeta)],
- [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]):
+ [diff(z_real, xi), diff(z_real, eta), diff(z_real, zeta)]]):
detJ := det (J):
measure := simplify ( int ( int ( int (detJ, xi=0..1), eta=0..1), zeta=0..1)):
For the case, someone is willing to rerun the maple script, he/she
should use following ordering of shape functions:
-
+
tphi[0] := (1-xi)*(1-eta)*(1-zeta):
tphi[1] := xi*(1-eta)*(1-zeta):
tphi[2] := (1-xi)* eta*(1-zeta):
accessor.vertex(6)(2) };
double s1, s2, s3, s4, s5, s6, s7, s8;
-
+
s1 = 1.0/6.0;
s8 = -x[2]*x[2]*y[0]*z[3]-2.0*z[6]*x[7]*x[7]*y[4]-z[5]*x[7]*x[7]*y[4]-z
[6]*x[7]*x[7]*y[5]+2.0*y[6]*x[7]*x[7]*z[4]-z[5]*x[6]*x[6]*y[4]+x[6]*x[6]*y[4]*z
\int_K 1 dx dy = \int_{\hat K} |det J| d(xi) d(eta)
# x and y are arrays holding the x- and y-values of the four vertices
- # of this cell in real space.
+ # of this cell in real space.
x := array(0..3);
y := array(0..3);
tphi[0] := (1-xi)*(1-eta):
#endif
-#if deal_II_dimension == 3
+#if deal_II_dimension == 3
double
measure (const TriaAccessor<3, 3, 3> &accessor)
}
#endif
-
+
template <int structdim, int dim, int spacedim>
double
template <>
double TriaAccessor<1,1,1>::extent_in_direction(const unsigned int axis) const
{
- Assert (axis == 0, ExcIndexRange (axis, 0, 1));
+ Assert (axis == 0, ExcIndexRange (axis, 0, 1));
return this->diameter();
}
-#elsif deal_II_dimension == 2
+#elif deal_II_dimension == 2
template <>
double TriaAccessor<2,2,2>::extent_in_direction(const unsigned int axis) const
{0,1}};/// Lines along y-axis
Assert (axis < 2, ExcIndexRange (axis, 0, 2));
-
+
return std::max(this->line(lines[axis][0])->diameter(),
this->line(lines[axis][1])->diameter());
}
-#elsif deal_II_dimension == 3
+#elif deal_II_dimension == 3
template <>
double TriaAccessor<3,3,3>::extent_in_direction(const unsigned int axis) const
const unsigned int lines[3][4] = {{2,3,6,7}, /// Lines along x-axis, see GeometryInfo
{0,1,4,5}, /// Lines along y-axis
{8,9,10,11}}; /// Lines along z-axis
-
+
Assert (axis < 3, ExcIndexRange (axis, 0, 3));
double lengths[4] = { this->line(lines[axis][0])->diameter(),
{
// original implementation by Joerg
// Weimar
-
+
// we first eliminate points based
// on the maximum and minumum of
// the corner coordinates, then
if (dim==1)
return GeometryInfo<dim>::opposite_face[neighbor];
-
+
const TriaIterator<CellAccessor<dim, spacedim> > neighbor_cell = this->neighbor(neighbor);
-
+
// usually, on regular patches of
// the grid, this cell is just on
// the opposite side of the
// GeometryInfo and try it
const unsigned int this_face_index=face_index(neighbor);
-
+
const unsigned int neighbor_guess
= GeometryInfo<dim>::opposite_face[neighbor];
-
+
if (neighbor_cell->face_index (neighbor_guess) == this_face_index)
return neighbor_guess;
else
// on a coarser level
Assert (neighbor_is_coarser(neighbor),
TriaAccessorExceptions::ExcNeighborIsNotCoarser());
-
+
switch (dim)
{
case 2:
{
const int this_face_index=face_index(neighbor);
const TriaIterator<CellAccessor<2,spacedim> > neighbor_cell = this->neighbor(neighbor);
-
+
// usually, on regular patches of
// the grid, this cell is just on
// the opposite side of the
const TriaIterator<TriaAccessor<1, 2, spacedim> > face_guess
=neighbor_cell->face(face_no_guess);
-
+
if (face_guess->has_children())
for (unsigned int subface_no=0; subface_no<face_guess->n_children(); ++subface_no)
if (face_guess->child_index(subface_no)==this_face_index)
return std::make_pair (face_no, subface_no);
}
}
-
+
// we should never get here,
// since then we did not find
// our way back...
const int this_face_index=face_index(neighbor);
const TriaIterator<CellAccessor<3, spacedim> >
neighbor_cell = this->neighbor(neighbor);
-
+
// usually, on regular patches of
// the grid, this cell is just on
// the opposite side of the
const TriaIterator<TriaAccessor<3-1, 3, spacedim> > face_guess
=neighbor_cell->face(face_no_guess);
-
+
if (face_guess->has_children())
for (unsigned int subface_no=0; subface_no<face_guess->n_children(); ++subface_no)
if (face_guess->child_index(subface_no)==this_face_index)
// a subface number for the current
// FaceRefineCase
return std::make_pair (face_no_guess, translate_subface_no(face_guess, subface_no, subsub_no));
-
-
+
+
// if the guess was false, then
// we need to loop over all faces
return std::make_pair (face_no, translate_subface_no(face, subface_no, subsub_no));
}
}
-
+
// we should never get here,
// since then we did not find
// our way back...
Assert (false, ExcInternalError());
return std::make_pair (numbers::invalid_unsigned_int,
- numbers::invalid_unsigned_int);
+ numbers::invalid_unsigned_int);
}
default:
Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
Assert (i<GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (i,0,GeometryInfo<dim>::faces_per_cell));
-
+
return (neighbor_index(i) == -1);
}
for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
if (this->line(l)->at_boundary())
return true;
-
+
return false;
}
}
ExcInternalError());
sub_neighbor = sub_neighbor->child(GeometryInfo<dim>::child_cell_on_face(
sub_neighbor->refinement_case(),neighbor_neighbor,0));
-
+
}
-
+
return sub_neighbor;
}
{
// this function returns the neighbor's
// child on a given face and
- // subface.
+ // subface.
// we have to consider one other aspect here:
// The face might be refined
//
// *---*---* *-------*
// | | | | 1 |
- // | 0 | 1 | or *-------*
+ // | 0 | 1 | or *-------*
// | | | | 0 |
// *---*---* *-------*
//
const unsigned int total_children=mother_face->number_of_children();
Assert (subface<total_children,ExcIndexRange(subface,0,total_children));
Assert (total_children<=GeometryInfo<3>::max_children_per_face, ExcInternalError());
-
+
unsigned int neighbor_neighbor;
TriaIterator<CellAccessor<3,spacedim> > neighbor_child;
const TriaIterator<CellAccessor<3,spacedim> > neighbor
= this->neighbor(face);
-
+
const RefinementCase<2> mother_face_ref_case
= mother_face->refinement_case();
if (mother_face_ref_case==RefinementCase<2>::cut_xy) // total_children==4
// face and the neighbor
neighbor_neighbor
= this->neighbor_of_neighbor (face);
-
+
// now use the info provided by GeometryInfo
// to extract the neighbors child number
const unsigned int neighbor_child_index
neighbor->face_flip(neighbor_neighbor),
neighbor->face_rotation(neighbor_neighbor));
neighbor_child = neighbor->child(neighbor_child_index);
-
+
// make sure that the neighbor child cell we
// have found shares the desired subface.
Assert((this->face(face)->child(subface) ==
unsigned int neighbor_child_index;
if (total_children==2)
first_child_to_find=subface;
- else
+ else
{
first_child_to_find=subface/2;
if (total_children==3 &&
std::pair<unsigned int, unsigned int> indices=neighbor_of_coarser_neighbor(face);
neighbor_neighbor=indices.first;
-
+
// we have to translate our
// subface_index according to the
// RefineCase and subface index of
neighbor->face_rotation(neighbor_neighbor),
mother_face_ref_case);
}
-
+
neighbor_child=neighbor->child(neighbor_child_index);
// it might be, that the neighbor_child
// has children, which are not refined
#endif
return neighbor_child;
-
+
}
default:
// 1d or more than 3d
Assert (false, ExcNotImplemented());
- return TriaIterator<CellAccessor<dim,spacedim> >();
+ return TriaIterator<CellAccessor<dim,spacedim> >();
}
}
projects / dealii-svn.git / commitdiff