source-clean:
cd source ; make clean
TAGS:
- etags include/*/* source/*/*
+ etags include/*/*.h source/*/*.cc
.PHONY: all deal.II examples TAGS
*/
const unsigned int n_transform_functions;
+
/**
- * Number of components.
+ * Number of components and dimension of the image space.
*/
const unsigned int n_components;
*
* The function assumes that the fields
* already have the right number of
- * elements.
+ * elements. It has to be
+ * guaranteed, that fields that are
+ * not requested for update are not changed.
+ * This also means, that these
+ * fields have to be filled with
+ * the correct values beforehand.
*
* This function is more or less an
* interface to the #FEValues# class and
const vector<Point<dim-1> > &unit_points,
vector<Point<dim> > &normal_vectors) const;
+ /**
+ * Implementation of the corresponding function of #FiniteElement#.
+ */
+ virtual void fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+ const vector<Point<dim> > &unit_points,
+ vector<Tensor<2,dim> > &jacobians,
+ const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
+ vector<Point<dim> > &q_points,
+ const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
+ const Boundary<dim> &boundary) const;
+
/**
* Number of different base
* elements of this object.
* the #.h# file.
*/
void initialize();
+ /**
+ *Exception.
+ */
+ DeclException0(ExcElementTransformNotEqual);
};
base_elements(2),
component_to_base_table(n_components)
{
+ Assert(fe1.n_transform_functions == fe2.n_transform_functions,
+ ExcElementTransformNotEqual());
+
base_elements[0] = ElementPair(new FE1, n1);
base_elements[0].first -> subscribe ();
base_elements[1] = ElementPair(new FE2, n2);
template <int dim> class StraightBoundary;
class ConstraintMatrix;
class dVector;
+class VectorFunction;
/**
map<int,double> &boundary_values);
/**
+ * Compute the error of the finite element solution.
* Integrate the difference between
* a finite element function and
* the reference function, which
const NormType &norm,
const Boundary<dim> &boundary=StraightBoundary<dim>());
+ /**
+ * Compute the error for the solution of a system.
+ * See the other #integrate_difference#.
+ */
+ static void integrate_difference (const DoFHandler<dim> &dof,
+ const dVector &fe_function,
+ const VectorFunction &exact_solution,
+ dVector &difference,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe,
+ const NormType &norm,
+ const Boundary<dim> &boundary);
+
+
/**
* Exception
*/
};
-
template<int dim>
bool FiniteElementData<dim>::operator== (const FiniteElementData<dim> &f) const
{
const dFMatrix &,
const vector<vector<Tensor<1,1> > > &,
const Boundary<1> &boundary) const {
- Assert (jacobians.size() == unit_points.size(),
+ Assert ((!compute_jacobians) || (jacobians.size() == unit_points.size()),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
- Assert (q_points.size() == unit_points.size(),
+ Assert ((!compute_jacobians_grad) || (jacobians_grad.size() == unit_points.size()),
+ ExcWrongFieldDimension(jacobians_grad.size(), unit_points.size()));
+ Assert ((!compute_q_points) || (q_points.size() == unit_points.size()),
ExcWrongFieldDimension(q_points.size(), unit_points.size()));
- Assert (support_points.size() == total_dofs,
+ Assert ((!compute_support_points) || (support_points.size() == total_dofs),
ExcWrongFieldDimension(support_points.size(), total_dofs));
vector<Tensor<3,dim> > &jacobians_grad,
const bool compute_jacobians_grad,
vector<Point<dim> > &support_points,
- const bool,
+ const bool compute_support_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
const dFMatrix &shape_values_transform,
// we need the support points in any
// way, wanted or not by the user
- get_support_points (cell, boundary, support_points);
+ if (compute_support_points)
+ get_support_points (cell, boundary, support_points);
if (compute_q_points)
{
const bool compute_q_points,
const dFMatrix &shape_values_transform,
const vector<vector<Tensor<1,dim> > > &/*shape_grad_transform*/,
- const Boundary<dim> &boundary) const {
- Assert (jacobians.size() == unit_points.size(),
+ const Boundary<dim> &boundary) const
+{
+ Assert ((!compute_jacobians) || (jacobians.size() == unit_points.size()),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
- Assert (q_points.size() == unit_points.size(),
+ Assert ((!compute_jacobians_grad) || (jacobians_grad.size() == unit_points.size()),
+ ExcWrongFieldDimension(jacobians_grad.size(), unit_points.size()));
+ Assert ((!compute_q_points) || (q_points.size() == unit_points.size()),
ExcWrongFieldDimension(q_points.size(), unit_points.size()));
- Assert (support_points.size() == total_dofs,
+ Assert ((!compute_support_points) || (support_points.size() == total_dofs),
ExcWrongFieldDimension(support_points.size(), total_dofs));
/* $Id$ */
-/* Copyright W. Bangerth, University of Heidelberg, 1990 */
+/* Copyright W. Bangerth, G. Kanschat University of Heidelberg, 1990 */
#include <fe/fe_system.h>
{
return FiniteElementData<1> (fe1.dofs_per_vertex * N1 + fe2.dofs_per_vertex * N2 ,
fe1.dofs_per_line * N1 + fe2.dofs_per_line * N2 ,
- fe1.n_transform_functions + fe2.n_transform_functions ,
+ fe1.n_transform_functions,
fe1.n_components * N1 + fe2.n_components * N2 );
};
return FiniteElementData<2> (fe1.dofs_per_vertex * N1 + fe2.dofs_per_vertex * N2 ,
fe1.dofs_per_line * N1 + fe2.dofs_per_line * N2 ,
fe1.dofs_per_quad * N1 + fe2.dofs_per_quad * N2 ,
- fe1.n_transform_functions + fe2.n_transform_functions ,
+ fe1.n_transform_functions,
fe1.n_components * N1 + fe2.n_components * N2 );
};
template <int dim>
-double FESystem<dim>::shape_value_transform (const unsigned int /*i*/,
- const Point<dim> &/*p*/) const
+double FESystem<dim>::shape_value_transform (const unsigned int i,
+ const Point<dim> &p) const
{
- Assert(false, ExcNotImplemented());
- return 0.;
+ return base_elements[0].first->shape_value_transform(i,p);
};
template <int dim>
-Tensor<1,dim> FESystem<dim>::shape_grad_transform (const unsigned int /*i*/,
- const Point<dim> &/*p*/) const
+Tensor<1,dim> FESystem<dim>::shape_grad_transform (const unsigned int i,
+ const Point<dim> &p) const
{
- Assert(false, ExcNotImplemented());
- return Tensor<1,dim>();
-
-// return base_element->shape_grad_transform (i, p);
+ return base_elements[0].first->shape_grad_transform (i, p);
};
template <int dim>
-void FESystem<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &/*face*/,
- const Boundary<dim> &/*boundary*/,
- const vector<Point<dim-1> > &/*unit_points*/,
- vector<double> &/*face_jacobi_determinants*/) const
+void FESystem<dim>::get_face_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const Boundary<dim> &boundary,
+ const vector<Point<dim-1> > &unit_points,
+ vector<double> &face_jacobi_determinants) const
{
- Assert(false, ExcNotImplemented());
-
-// base_element->get_face_jacobians (face, boundary, unit_points, face_jacobi_determinants);
+ base_elements[0].first->get_face_jacobians (face, boundary, unit_points,
+ face_jacobi_determinants);
};
template <int dim>
-void FESystem<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &/*face*/,
- const unsigned int /*subface_no*/,
- const vector<Point<dim-1> > &/*unit_points*/,
- vector<double> &/*face_jacobi_determinants*/) const
+void FESystem<dim>::get_subface_jacobians (const DoFHandler<dim>::face_iterator &face,
+ const unsigned int subface_no,
+ const vector<Point<dim-1> > &unit_points,
+ vector<double> &face_jacobi_determinants) const
{
- Assert(false, ExcNotImplemented());
-
-// base_element->get_subface_jacobians (face, subface_no, unit_points, face_jacobi_determinants);
+ base_elements[0].first->get_subface_jacobians (face, subface_no, unit_points,
+ face_jacobi_determinants);
};
template <int dim>
-void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &/*cell*/,
- const unsigned int /*face_no*/,
- const Boundary<dim> &/*boundary*/,
- const vector<Point<dim-1> > &/*unit_points*/,
- vector<Point<dim> > &/*normal_vectors*/) const
+void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Boundary<dim> &boundary,
+ const vector<Point<dim-1> > &unit_points,
+ vector<Point<dim> > &normal_vectors) const
{
- Assert(false, ExcNotImplemented());
-
-// base_element->get_normal_vectors (cell, face_no, boundary, unit_points, normal_vectors);
+ base_elements[0].first->get_normal_vectors (cell, face_no, boundary, unit_points,
+ normal_vectors);
};
template <int dim>
-void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &/*cell*/,
- const unsigned int /*face_no*/,
- const unsigned int /*subface_no*/,
- const vector<Point<dim-1> > &/*unit_points*/,
- vector<Point<dim> > &/*normal_vectors*/) const
+void FESystem<dim>::get_normal_vectors (const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const vector<Point<dim-1> > &unit_points,
+ vector<Point<dim> > &normal_vectors) const
{
- Assert(false, ExcNotImplemented());
-
-// base_element->get_normal_vectors (cell, face_no, subface_no, unit_points, normal_vectors);
+ base_elements[0].first->get_normal_vectors (cell, face_no, subface_no, unit_points,
+ normal_vectors);
};
+template <int dim>
+void
+FESystem<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
+ const vector<Point<dim> > &unit_points,
+ vector<Tensor<2,dim> > &jacobians,
+ const bool compute_jacobians,
+ vector<Tensor<3,dim> > &jacobians_grad,
+ const bool compute_jacobians_grad,
+ vector<Point<dim> > &support_points,
+ const bool compute_support_points,
+ vector<Point<dim> > &q_points,
+ const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Tensor<1,dim> > > &shape_grad_transform,
+ const Boundary<dim> &boundary) const
+{
+ vector<Point<dim> > supp(base_elements[0].first->total_dofs);
+ base_elements[0].first->fill_fe_values (cell, unit_points, jacobians, compute_jacobians,
+ jacobians_grad, compute_jacobians_grad,
+ support_points, compute_support_points,
+ q_points, compute_q_points,
+ shape_values_transform, shape_grad_transform, boundary);
+
+ if (compute_support_points)
+ {
+ unsigned component = 0;
+
+ for (unsigned m=0 ; m < element_multiplicity(0) ; ++ m)
+ {
+ for (unsigned i=0 ; i < base_element(0).total_dofs ; ++i)
+ support_points[component_to_system_index(component,i)] = supp[i];
+ ++component;
+ }
+ for (unsigned base=1 ; base < n_base_elements() ; ++base)
+ {
+ supp.resize(base_elements[base].first->total_dofs);
+ base_elements[base].first->fill_fe_values (cell, unit_points, jacobians, false,
+ jacobians_grad, false,
+ supp, true,
+ q_points, false,
+ shape_values_transform, shape_grad_transform, boundary);
+
+ for (unsigned m=0 ; m < element_multiplicity(base) ; ++ m)
+ {
+ for (unsigned i=0 ; i < base_element(base).total_dofs ; ++i)
+ support_points[component_to_system_index(component,i)] = supp[i];
+ ++component;
+ }
+ }
+ }
+}
+
// explicit instantiations
template class FESystem<deal_II_dimension>;
};
+template <int dim>
+void FEValuesBase<dim>::get_function_values (const dVector &fe_function,
+ vector< vector<double> > &values) const
+{
+ Assert (fe->n_components == values.size(),
+ ExcWrongNoOfComponents());
+ Assert (selected_dataset<shape_values.size(),
+ ExcInvalidIndex (selected_dataset, shape_values.size()));
+ for (unsigned i=0;i<values.size();++i)
+ Assert (values[i].size() == n_quadrature_points,
+ ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+ // get function values of dofs
+ // on this cell
+ dVector dof_values (total_dofs);
+ present_cell->get_dof_values (fe_function, dof_values);
+
+ // initialize with zero
+ for (unsigned i=0;i<values.size();++i)
+ fill_n (values[i].begin(), n_quadrature_points, 0);
+
+ // add up contributions of trial
+ // functions
+ for (unsigned int point=0; point<n_quadrature_points; ++point)
+ for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
+ values[fe->system_to_component_index(shape_func).first][point]
+ += (dof_values(shape_func) * shape_values[selected_dataset](shape_func, point));
+};
+
+
template <int dim>
const Tensor<1,dim> &
FEValuesBase<dim>::shape_grad (const unsigned int i,
- const unsigned int j) const {
+ const unsigned int j) const
+{
Assert (i<shape_gradients.size(),
ExcInvalidIndex (i, shape_gradients.size()));
Assert (j<shape_gradients[i].size(),
FEValues<dim> fe(dofs->get_fe(), points, UpdateFlags(update_q_points));
const StraightBoundary<dim> boundary;
- vector<double> *values = new vector<double> [data.size()];
+ vector< vector <vector<double> > >
+ values (data.size(), vector< vector<double> >(dofs->get_fe().n_components, vector<double>(points.n_quadrature_points)));
for (cell=dofs->begin_active(); cell!=endc; ++cell)
{
for (unsigned i=0; i<data.size(); ++i)
{
- values[i].resize(points.n_quadrature_points);
+// values[i].resize(points.n_quadrature_points);
fe.get_function_values(*data[i].data, values[i]);
}
Point<dim> pt = fe.quadrature_point(supp_pt);
out << pt << " ";
for (unsigned int i=0; i!=data.size(); ++i)
- out << values[i][supp_pt]
- << ' ';
+ for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
+ out << values[i][j][supp_pt]
+ << ' ';
out << endl;
};
out << pt << " ";
for (unsigned int i=0; i!=data.size(); ++i)
- out << values[i][supp_pt]
- << ' ';
+ for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
+ out << values[i][j][supp_pt]
+ << ' ';
out << endl;
}
out << endl;
out << pt << " ";
for (unsigned int i=0; i!=data.size(); ++i)
- out << values[i][supp_pt]
- << ' ';
+ for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
+ out << values[i][j][supp_pt]
+ << ' ';
out << endl;
}
out << endl;
}
out << endl;
}
- delete [] values;
AssertThrow (out, ExcIO());
}
};
};
+template <int dim>
+void VectorTools<dim>::integrate_difference (const DoFHandler<dim> &dof,
+ const dVector &fe_function,
+ const TensorFunction<1,dim> &exact_solution,
+ dVector &difference,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe,
+ const NormType &norm,
+ const Boundary<dim> &boundary) {
+ Assert (fe == dof.get_fe(), ExcInvalidFE());
+
+ difference.reinit (dof.get_tria().n_active_cells());
+
+ UpdateFlags update_flags = UpdateFlags (update_q_points |
+ update_JxW_values);
+ if ((norm==H1_seminorm) || (norm==H1_norm))
+ update_flags = UpdateFlags (update_flags | update_gradients);
+ FEValues<dim> fe_values(fe, q, update_flags);
+
+ // loop over all cells
+ DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+ endc = dof.end();
+ for (unsigned int index=0; cell != endc; ++cell, ++index)
+ {
+ double diff=0;
+ // initialize for this cell
+ fe_values.reinit (cell, boundary);
+
+ switch (norm)
+ {
+ case mean:
+ case L1_norm:
+ case L2_norm:
+ case Linfty_norm:
+ case H1_norm:
+ {
+ // we need the finite element
+ // function \psi at the different
+ // integration points. Compute
+ // it like this:
+ // \psi(x_j)=\sum_i v_i \phi_i(x_j)
+ // with v_i the nodal values of the
+ // fe_function and \phi_i(x_j) the
+ // matrix of the trial function
+ // values at the integration point
+ // x_j. Then the vector
+ // of the \psi(x_j) is v*Phi with
+ // v being the vector of nodal
+ // values on this cell and Phi
+ // the matrix.
+ //
+ // we then need the difference:
+ // reference_function(x_j)-\psi_j
+ // and assign that to the vector
+ // \psi.
+ const unsigned int n_q_points = q.n_quadrature_points;
+ vector<double> psi (n_q_points);
+
+ // in praxi: first compute
+ // exact fe_function vector
+ exact_solution.value_list (fe_values.get_quadrature_points(),
+ psi);
+ // then subtract finite element
+ // fe_function
+ if (true)
+ {
+ vector< vector<double> > function_values (fe.n_components, vector<double>(n_q_points, 0));
+ fe_values.get_function_values (fe_function, function_values);
+
+ transform (psi.begin(), psi.end(),
+ function_values.begin(),
+ psi.begin(),
+ minus<double>());
+ };
+
+ // for L1_norm and Linfty_norm:
+ // take absolute
+ // value, for the L2_norm take
+ // square of psi
+ switch (norm)
+ {
+ case mean:
+ break;
+ case L1_norm:
+ case Linfty_norm:
+ transform (psi.begin(), psi.end(),
+ psi.begin(), ptr_fun(fabs));
+ break;
+ case L2_norm:
+ case H1_norm:
+ transform (psi.begin(), psi.end(),
+ psi.begin(), ptr_fun(sqr));
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ // ok, now we have the integrand,
+ // let's compute the integral,
+ // which is
+ // sum_j psi_j JxW_j
+ // (or |psi_j| or |psi_j|^2
+ switch (norm)
+ {
+ case mean:
+ case L1_norm:
+ diff = inner_product (psi.begin(), psi.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0);
+ break;
+ case L2_norm:
+ case H1_norm:
+ diff = sqrt(inner_product (psi.begin(), psi.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0));
+ break;
+ case Linfty_norm:
+ diff = *max_element (psi.begin(), psi.end());
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ // note: the H1_norm uses the result
+ // of the L2_norm and control goes
+ // over to the next case statement!
+ if (norm != H1_norm)
+ break;
+ };
+
+ case H1_seminorm:
+ {
+ // note: the computation of the
+ // H1_norm starts at the previous
+ // case statement, but continues
+ // here!
+
+ // for H1_norm: re-square L2_norm.
+ diff = sqr(diff);
+
+ // same procedure as above, but now
+ // psi is a vector of gradients
+ const unsigned int n_q_points = q.n_quadrature_points;
+ vector<Tensor<1,dim> > psi (n_q_points);
+
+ // in praxi: first compute
+ // exact fe_function vector
+ exact_solution.gradient_list (fe_values.get_quadrature_points(),
+ psi);
+
+ // then subtract finite element
+ // fe_function
+ if (true)
+ {
+ vector<Tensor<1,dim> > function_grads (n_q_points, Tensor<1,dim>());
+ fe_values.get_function_grads (fe_function, function_grads);
+
+ transform (psi.begin(), psi.end(),
+ function_grads.begin(),
+ psi.begin(),
+ minus<Tensor<1,dim> >());
+ };
+ // take square of integrand
+ vector<double> psi_square (psi.size(), 0.0);
+ for (unsigned int i=0; i<n_q_points; ++i)
+ psi_square[i] = sqr_point(psi[i]);
+
+ // add seminorm to L_2 norm or
+ // to zero
+ diff += inner_product (psi_square.begin(), psi_square.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0);
+ diff = sqrt(diff);
+
+ break;
+ };
+
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+
+ // append result of this cell
+ // to the end of the vector
+ difference(index) = diff;
+ };
+};
+
template VectorTools<deal_II_dimension>;