delete base_fe_datas[i];
base_fe_datas[i] = 0;
};
-
+
for (unsigned int i=0; i<base_fe_values_datas.size(); ++i)
if (base_fe_values_datas[i])
{
fe3, n3)),
base_elements(3)
{
- base_elements[0] = ElementPair(fe1.clone(), n1);
+ base_elements[0] = ElementPair(fe1.clone(), n1);
base_elements[0].first->subscribe (typeid(*this).name());
base_elements[1] = ElementPair(fe2.clone(), n2);
base_elements[1].first->subscribe (typeid(*this).name());
const Point<dim> &p) const
{
Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- Assert (this->is_primitive(i),
+ Assert (this->is_primitive(i),
(typename FiniteElement<dim,spacedim>::ExcShapeFunctionNotPrimitive(i)));
return (base_element(this->system_to_base_table[i].first.first)
// zero, then return right away...
if (this->nonzero_components[i][component] == false)
return 0;
-
+
// ...otherwise: first find out to
// which of the base elements this
// desired component belongs, and
Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
Assert (component < this->n_components(),
ExcIndexRange (component, 0, this->n_components()));
-
+
// if this value is supposed to be
// zero, then return right away...
if (this->nonzero_components[i][component] == false)
// element it is
const unsigned int base = this->component_to_base_index(component).first;
const unsigned int component_in_base = this->component_to_base_index(component).second;
-
+
// then get value from base
// element. note that that will
// throw an error should the
const Point<dim> &p) const
{
Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
- Assert (this->is_primitive(i),
+ Assert (this->is_primitive(i),
(typename FiniteElement<dim,spacedim>::ExcShapeFunctionNotPrimitive(i)));
return (base_element(this->system_to_base_table[i].first.first)
Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
Assert (component < this->n_components(),
ExcIndexRange (component, 0, this->n_components()));
-
+
// if this value is supposed to be
// zero, then return right away...
if (this->nonzero_components[i][component] == false)
// element it is
const unsigned int base = this->component_to_base_index(component).first;
const unsigned int component_in_base = this->component_to_base_index(component).second;
-
+
// then get value from base
// element. note that that will
// throw an error should the
// there are certain conditions
// that the two elements have to
// satisfy so that this can work.
- //
+ //
// condition 1: the other element
// must also be a system element
(dynamic_cast<const FESystem<dim,spacedim>*>(&x_source_fe) != 0),
typename FEL::
ExcInterpolationNotImplemented());
-
+
// ok, source is a system element,
// so we may be able to do the work
const FESystem<dim,spacedim> &source_fe
// ok, so let's try whether it
// works:
-
+
// first let's see whether all the
// basis elements actually generate
// their interpolation matrices. if
// elements don't have them
if (flags & update_hessians)
out |= update_hessians | update_covariant_transformation;
-
+
return out;
}
const Quadrature<dim> &quadrature) const
{
UpdateFlags flags = flags_;
-
+
InternalData* data = new InternalData(n_base_elements());
data->update_once = update_once (flags);
data->update_each = update_each (flags);
flags = data->update_once | data->update_each;
-
+
UpdateFlags sub_flags = flags;
// if second derivatives through
sub_flags = UpdateFlags (sub_flags ^ update_hessians);
data->initialize_2nd (this, mapping, quadrature);
}
-
+
// get data objects from each of
// the base elements and store them
typename FiniteElement<dim,spacedim>::InternalDataBase *base_fe_data =
dynamic_cast<typename FiniteElement<dim,spacedim>::InternalDataBase *>
(base_fe_data_base);
-
+
data->set_fe_data(base_no, base_fe_data);
// make sure that *we* compute
ExcInternalError());
Assert (!(base_fe_data->update_once & update_hessians),
ExcInternalError());
-
+
// The FEValuesData @p{data}
// given to the
// @p{fill_fe_values} function
typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
typename Mapping<dim,spacedim>::InternalDataBase &fedata,
FEValuesData<dim,spacedim> &data) const
-{
+{
const unsigned int n_q_points = quadrature.size();
-
+
// convert data object to internal
// data for this class. fails with
// an exception if that is not
// possible
Assert (dynamic_cast<InternalData*> (&fedata) != 0, ExcInternalError());
InternalData & fe_data = static_cast<InternalData&> (fedata);
-
+
// Either dim_1==dim
// (fill_fe_values) or dim_1==dim-1
// (fill_fe_(sub)face_values)
// base element
typename Mapping<dim,spacedim>::InternalDataBase &
base_fe_data = fe_data.get_fe_data(base_no);
-
+
// compute update flags ...
const UpdateFlags base_update_flags
= mapping_data.update_flags | base_fe_data.update_flags;
-
+
// Initialize the FEValuesDatas
// for the base elements.
FEValuesData<dim,spacedim> &base_data=fe_data.get_fe_values_data(base_no);
base_data.initialize (n_q_points, base_fe, base_update_flags);
}
}
-
+
// fill_fe_face_values needs
// argument Quadrature<dim-1>
// for both cases
// following workaround
const Quadrature<dim> *cell_quadrature = 0;
const Quadrature<dim-1> *face_quadrature = 0;
-
+
// static cast to the
// common base class of
// quadrature being either
if (face_no==invalid_face_number)
{
Assert(dim_1==dim, ExcDimensionMismatch(dim_1,dim));
- Assert (dynamic_cast<const Quadrature<dim> *>(quadrature_base_pointer) != 0,
+
+ // Some Mac OS X versions
+ // have a bug that makes
+ // the following assertion
+ // fail spuriously. In that
+ // case, just disable the
+ // Assert, as this is not
+ // something that depends
+ // on input parameters
+ // provided by users; the
+ // validity of the
+ // assertion is still
+ // checked on all systems
+ // that don't have this bug
+#ifndef DEAL_II_HAVE_DARWIN_DYNACAST_BUG
+ Assert (dynamic_cast<const Quadrature<dim> *>(quadrature_base_pointer) != 0,
ExcInternalError());
+#endif
cell_quadrature
= static_cast<const Quadrature<dim> *>(quadrature_base_pointer);
}
else
{
Assert(dim_1==dim-1, ExcDimensionMismatch(dim_1,dim-1));
- Assert (dynamic_cast<const Quadrature<dim-1> *>(quadrature_base_pointer) != 0,
+#ifndef DEAL_II_HAVE_DARWIN_DYNACAST_BUG
+ Assert (dynamic_cast<const Quadrature<dim-1> *>(quadrature_base_pointer) != 0,
ExcInternalError());
+#endif
face_quadrature
= static_cast<const Quadrature<dim-1> *>(quadrature_base_pointer);
}
// therefore use
// base_fe_data.update_flags.
if (face_no==invalid_face_number)
- base_fe.fill_fe_values(mapping, cell, *cell_quadrature, mapping_data,
+ base_fe.fill_fe_values(mapping, cell, *cell_quadrature, mapping_data,
base_fe_data, base_data, cell_similarity);
else if (sub_no==invalid_face_number)
base_fe.fill_fe_face_values(mapping, cell, face_no,
unsigned int in_index = 0;
for (unsigned int i=0; i<base_index; ++i)
in_index += base_fe.n_nonzero_components(i);
-
+
// then loop over the
// number of components
// to be copied
for (unsigned int q=0; q<n_q_points; ++q)
data.shape_values[out_index+s][q] =
base_data.shape_values(in_index+s,q);
-
+
if (base_flags & update_gradients)
for (unsigned int q=0; q<n_q_points; ++q)
data.shape_gradients[out_index+s][q]=
// base element
typename Mapping<dim,spacedim>::InternalDataBase &base_fe_data=
fe_data.get_fe_data(base_no);
-
+
// compute update flags ...
UpdateFlags base_flags_each(
dim_1==dim ?
mapping_data.update_each | base_fe_data.update_each :
mapping_data.update_flags | base_fe_data.update_flags);
-
+
if (base_flags_each==update_default)
fe_data.delete_fe_values_data(base_no);
}
}
}
-
+
if (fe_data.compute_hessians && cell_similarity != CellSimilarity::translation)
{
unsigned int offset = 0;
n_q_points,
cell->subface_case(face_no));
}
-
+
this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
}
}
// by FiniteElement
this->system_to_component_table.resize(this->dofs_per_cell);
this->face_system_to_component_table.resize(this->dofs_per_face);
-
+
unsigned total_index = 0;
for (unsigned int base=0; base < n_base_elements(); ++base)
for (unsigned int m = 0; m < element_multiplicity(base); ++m)
const std::pair<unsigned int, unsigned int>
non_primitive_index (numbers::invalid_unsigned_int,
numbers::invalid_unsigned_int);
-
+
// First enumerate vertex indices,
// where we first enumerate all
// indices on the first vertex in
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (base_element(base).dofs_per_vertex*vertex_number +
+ = (base_element(base).dofs_per_vertex*vertex_number +
local_index);
this->system_to_base_table[total_index]
this->system_to_component_table[total_index] = non_primitive_index;
}
}
-
+
// 2. Lines
if (GeometryInfo<dim>::lines_per_cell > 0)
for (unsigned int line_number= 0;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (base_element(base).dofs_per_line*line_number +
+ = (base_element(base).dofs_per_line*line_number +
local_index +
base_element(base).first_line_index);
-
+
this->system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), index_in_base);
-
+
if (base_element(base).is_primitive(index_in_base))
{
const unsigned int comp_in_base
this->system_to_component_table[total_index] = non_primitive_index;
}
}
-
+
// 3. Quads
if (GeometryInfo<dim>::quads_per_cell > 0)
for (unsigned int quad_number= 0;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (base_element(base).dofs_per_quad*quad_number +
+ = (base_element(base).dofs_per_quad*quad_number +
local_index +
base_element(base).first_quad_index);
-
+
this->system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), index_in_base);
-
+
if (base_element(base).is_primitive(index_in_base))
{
const unsigned int comp_in_base
this->system_to_component_table[total_index] = non_primitive_index;
}
}
-
+
// 4. Hexes
if (GeometryInfo<dim>::hexes_per_cell > 0)
for (unsigned int hex_number= 0;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (base_element(base).dofs_per_hex*hex_number +
+ = (base_element(base).dofs_per_hex*hex_number +
local_index +
base_element(base).first_hex_index);
-
+
this->system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), index_in_base);
-
+
if (base_element(base).is_primitive(index_in_base))
{
const unsigned int comp_in_base
const std::pair<unsigned int, unsigned int>
non_primitive_index (numbers::invalid_unsigned_int,
numbers::invalid_unsigned_int);
-
+
// 1. Vertices
unsigned int total_index = 0;
for (unsigned int vertex_number=0;
// on the same type of
// object):
const unsigned int index_in_base
- = (base_element(base).dofs_per_vertex*vertex_number +
+ = (base_element(base).dofs_per_vertex*vertex_number +
local_index);
-
+
const unsigned int face_index_in_base
- = (base_element(base).dofs_per_vertex*vertex_number +
+ = (base_element(base).dofs_per_vertex*vertex_number +
local_index);
this->face_system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), face_index_in_base);
-
+
if (base_element(base).is_primitive(index_in_base))
{
const unsigned int comp_in_base
this->face_system_to_component_table[total_index] = non_primitive_index;
}
}
-
+
// 2. Lines
if (GeometryInfo<dim>::lines_per_face > 0)
for (unsigned line_number= 0;
// alike for this
// type of object
const unsigned int index_in_base
- = (base_element(base).dofs_per_line*line_number +
+ = (base_element(base).dofs_per_line*line_number +
local_index +
base_element(base).first_line_index);
-
+
const unsigned int face_index_in_base
= (base_element(base).first_face_line_index +
- base_element(base).dofs_per_line * line_number +
+ base_element(base).dofs_per_line * line_number +
local_index);
this->face_system_to_base_table[total_index]
this->face_system_to_component_table[total_index] = non_primitive_index;
}
}
-
+
// 3. Quads
if (GeometryInfo<dim>::quads_per_face > 0)
for (unsigned quad_number= 0;
// alike for this
// type of object
const unsigned int index_in_base
- = (base_element(base).dofs_per_quad*quad_number +
+ = (base_element(base).dofs_per_quad*quad_number +
local_index +
base_element(base).first_quad_index);
-
+
const unsigned int face_index_in_base
= (base_element(base).first_face_quad_index +
- base_element(base).dofs_per_quad * quad_number +
+ base_element(base).dofs_per_quad * quad_number +
local_index);
-
+
this->face_system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), face_index_in_base);
template <int dim, int spacedim>
-void FESystem<dim,spacedim>::build_interface_constraints ()
+void FESystem<dim,spacedim>::build_interface_constraints ()
{
// check whether all base elements
// implement their interface
for (unsigned int base=0; base<n_base_elements(); ++base)
if (base_element(base).constraints_are_implemented() == false)
return;
-
+
this->interface_constraints.
TableBase<2,double>::reinit (this->interface_constraints_size());
-
+
// the layout of the constraints
// matrix is described in the
// FiniteElement class. you may
// same way as above,
// although a little
// more complicated...
-
+
// the indices
// m=0..5*d_v-1 are
// from the center and
base_element(m_index.first.first).dofs_per_quad*sub_quad +
tmp2;
};
-
+
break;
};
-
+
default:
Assert (false, ExcNotImplemented());
};
if (base_element(i).prolongation[ref_case-1][0].n() == 0)
do_prolongation = false;
}
-
+
// if we did not encounter void
// matrices, initialize the
// respective matrix sizes
for (unsigned int i=0;i<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));++i)
this->prolongation[ref_case-1][i].reinit(this->dofs_per_cell,
this->dofs_per_cell);
-
+
// distribute the matrices of the
// base finite elements to the
// matrices of this object. for
= (base_element(base)
.get_restriction_matrix(child, RefinementCase<dim>(ref_case))(
base_index_i,base_index_j));
-
+
if (do_prolongation)
this->prolongation[ref_case-1][child] (i,j)
= (base_element(base)
dpo.push_back(fe_data.dofs_per_line * N);
if (dim>1) dpo.push_back(fe_data.dofs_per_quad * N);
if (dim>2) dpo.push_back(fe_data.dofs_per_hex * N);
-
+
return FiniteElementData<dim> (dpo,
fe_data.n_components() * N,
fe_data.tensor_degree(),
this->unit_face_support_points.resize(0);
return;
};
-
+
// generate unit face support points
// from unit support points of sub
Assert (index_in_base < base_element(base_i).unit_face_support_points.size(),
ExcInternalError());
-
+
this->unit_face_support_points[i]
= base_element(base_i).unit_face_support_points[index_in_base];
};
// clear matrix, since we will not get to
// set all elements
interpolation_matrix = 0;
-
+
// loop over all the base elements of
// this and the other element, counting
// their multiplicities
multiplicity_other = 0;
FullMatrix<double> base_to_base_interpolation;
-
+
while (true)
{
const FiniteElement<dim,spacedim>
base.dofs_per_face);
base.get_face_interpolation_matrix (base_other,
base_to_base_interpolation);
-
+
// now translate entries. we'd like to
// have something like
// face_base_to_system_index, but that
interpolation_matrix(j, i)
= base_to_base_interpolation(fe_other_system->face_system_to_base_index(j).second,
this->face_system_to_base_index(i).second);
-
+
// advance to the next base element
// for this and the other
// fe_system; see if we can simply
// clear matrix, since we will not get to
// set all elements
interpolation_matrix = 0;
-
+
// loop over all the base elements of
// this and the other element, counting
// their multiplicities
multiplicity_other = 0;
FullMatrix<double> base_to_base_interpolation;
-
+
while (true)
{
const FiniteElement<dim,spacedim>
base.get_subface_interpolation_matrix (base_other,
subface,
base_to_base_interpolation);
-
+
// now translate entries. we'd like to
// have something like
// face_base_to_system_index, but that
interpolation_matrix(j, i)
= base_to_base_interpolation(fe_other_system->face_system_to_base_index(j).second,
this->face_system_to_base_index(i).second);
-
+
// advance to the next base element
// for this and the other
// fe_system; see if we can simply
dof_offset_other = 0;
std::vector<std::pair<unsigned int, unsigned int> > identities;
-
+
while (true)
{
const FiniteElement<dim,spacedim>
break;
default:
Assert (false, ExcNotImplemented());
- }
+ }
for (unsigned int i=0; i<base_identities.size(); ++i)
identities.push_back
// already taken care of
dof_offset += base.template n_dofs_per_object<structdim>();
dof_offset_other += base_other.template n_dofs_per_object<structdim>();
-
+
// advance to the next base element
// for this and the other
// fe_system; see if we can simply
Assert (domination !=
FiniteElementDomination::neither_element_dominates,
ExcInternalError());
-
+
return domination;
}
-
+
Assert (false, ExcNotImplemented());
return FiniteElementDomination::neither_element_dominates;
}
template <int dim, int spacedim>
std::vector<bool>
FESystem<dim,spacedim>::compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> &fe,
- const unsigned int n_elements)
+ const unsigned int n_elements)
{
std::vector<const FiniteElement<dim,spacedim>*> fe_list;
std::vector<unsigned int> multiplicities;
fe_list.push_back (&fe);
multiplicities.push_back (n_elements);
-
+
return compute_restriction_is_additive_flags (fe_list, multiplicities);
}
FESystem<dim,spacedim>::compute_restriction_is_additive_flags (const FiniteElement<dim,spacedim> &fe1,
const unsigned int N1,
const FiniteElement<dim,spacedim> &fe2,
- const unsigned int N2)
+ const unsigned int N2)
{
std::vector<const FiniteElement<dim,spacedim>*> fe_list;
std::vector<unsigned int> multiplicities;
fe_list.push_back (&fe2);
multiplicities.push_back (N2);
-
+
return compute_restriction_is_additive_flags (fe_list, multiplicities);
}
const FiniteElement<dim,spacedim> &fe2,
const unsigned int N2,
const FiniteElement<dim,spacedim> &fe3,
- const unsigned int N3)
+ const unsigned int N3)
{
std::vector<const FiniteElement<dim,spacedim>*> fe_list;
std::vector<unsigned int> multiplicities;
fe_list.push_back (&fe3);
multiplicities.push_back (N3);
-
+
return compute_restriction_is_additive_flags (fe_list, multiplicities);
}
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (fes[base]->dofs_per_vertex*vertex_number +
+ = (fes[base]->dofs_per_vertex*vertex_number +
local_index);
Assert (index_in_base < fes[base]->dofs_per_cell,
retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
}
}
-
+
// 2. Lines
if (GeometryInfo<dim>::lines_per_cell > 0)
for (unsigned int line_number= 0;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (fes[base]->dofs_per_line*line_number +
+ = (fes[base]->dofs_per_line*line_number +
local_index +
fes[base]->first_line_index);
retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
}
}
-
+
// 3. Quads
if (GeometryInfo<dim>::quads_per_cell > 0)
for (unsigned int quad_number= 0;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (fes[base]->dofs_per_quad*quad_number +
+ = (fes[base]->dofs_per_quad*quad_number +
local_index +
fes[base]->first_quad_index);
retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
}
}
-
+
// 4. Hexes
if (GeometryInfo<dim>::hexes_per_cell > 0)
for (unsigned int hex_number= 0;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (fes[base]->dofs_per_hex*hex_number +
+ = (fes[base]->dofs_per_hex*hex_number +
local_index +
fes[base]->first_hex_index);
}
Assert (total_index == n_shape_functions, ExcInternalError());
-
+
return retval;
}
fe_list.push_back (&fe1);
multiplicities.push_back (N1);
-
+
return compute_nonzero_components (fe_list, multiplicities);
}
fe_list.push_back (&fe2);
multiplicities.push_back (N2);
-
+
return compute_nonzero_components (fe_list, multiplicities);
}
fe_list.push_back (&fe3);
multiplicities.push_back (N3);
-
+
return compute_nonzero_components (fe_list, multiplicities);
}
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (fes[base]->dofs_per_vertex*vertex_number +
+ = (fes[base]->dofs_per_vertex*vertex_number +
local_index);
Assert (comp_start+fes[base]->n_components() <=
};
}
}
-
+
// 2. Lines
if (GeometryInfo<dim>::lines_per_cell > 0)
for (unsigned int line_number= 0;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (fes[base]->dofs_per_line*line_number +
+ = (fes[base]->dofs_per_line*line_number +
local_index +
fes[base]->first_line_index);
};
}
}
-
+
// 3. Quads
if (GeometryInfo<dim>::quads_per_cell > 0)
for (unsigned int quad_number= 0;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (fes[base]->dofs_per_quad*quad_number +
+ = (fes[base]->dofs_per_quad*quad_number +
local_index +
fes[base]->first_quad_index);
};
}
}
-
+
// 4. Hexes
if (GeometryInfo<dim>::hexes_per_cell > 0)
for (unsigned int hex_number= 0;
++local_index, ++total_index)
{
const unsigned int index_in_base
- = (fes[base]->dofs_per_hex*hex_number +
+ = (fes[base]->dofs_per_hex*hex_number +
local_index +
fes[base]->first_hex_index);
}
Assert (total_index == n_shape_functions, ExcInternalError());
-
+
return retval;
}
const FiniteElement<dim,spacedim> &
FESystem<dim,spacedim>::base_element (const unsigned int index) const
{
- Assert (index < base_elements.size(),
+ Assert (index < base_elements.size(),
ExcIndexRange(index, 0, base_elements.size()));
return *base_elements[index].first;
}
unsigned int
FESystem<dim,spacedim>::element_multiplicity (const unsigned int index) const
{
- Assert (index < base_elements.size(),
+ Assert (index < base_elements.size(),
ExcIndexRange(index, 0, base_elements.size()));
return base_elements[index].second;
}
template <int dim, int spacedim>
unsigned int
-FESystem<dim,spacedim>::memory_consumption () const
+FESystem<dim,spacedim>::memory_consumption () const
{
// neglect size of data stored in
// @p{base_elements} due to some