// children.
fine.reinit(tria.begin_active());
FullMatrix<number> A(nq*nd, n);
- for (unsigned int d=0;d<nd;++d)
- for (unsigned int k=0;k<nq;++k)
- for (unsigned int j=0;j<n;++j)
+ for (unsigned int j=0;j<n;++j)
+ for (unsigned int d=0;d<nd;++d)
+ for (unsigned int k=0;k<nq;++k)
A(k*nd+d,j) = fine.shape_value_component(j,k,d);
Householder<double> H(A);
coarse.reinit(tria.begin(0));
FullMatrix<double> &this_matrix = matrices[ref_case-1][cell_number];
+ v_coarse = 0;
// Compute this once for each
// coarse grid basis function
// function values of the
// coarse grid function in
// each quadrature point.
- for (unsigned int d=0;d<nd;++d)
- for (unsigned int k=0;k<nq;++k)
- v_coarse(k*nd+d) = coarse.shape_value_component(i,k,d);
+ if (fe.is_primitive())
+ {
+ const unsigned int d = fe.system_to_component_index(i).first;
+ const double * phi_i = &coarse.shape_value (i,0);
+ for (unsigned int k=0;k<nq;++k)
+ v_coarse(k*nd+d) = phi_i[k];
+ }
+ else
+ for (unsigned int d=0;d<nd;++d)
+ for (unsigned int k=0;k<nq;++k)
+ v_coarse(k*nd+d) = coarse.shape_value_component(i,k,d);
// solve the least squares
// problem.
// children.
fine.reinit(tria.begin_active());
FullMatrix<number> A(nq*nd, n);
- for (unsigned int d=0;d<nd;++d)
- for (unsigned int k=0;k<nq;++k)
- for (unsigned int j=0;j<n;++j)
+ for (unsigned int j=0;j<n;++j)
+ for (unsigned int d=0;d<nd;++d)
+ for (unsigned int k=0;k<nq;++k)
A(k*nd+d,j) = fine.shape_value_component(face_f_dofs[j],k,d);
Householder<double> H(A);
// each quadrature point.
for (unsigned int d=0;d<nd;++d)
for (unsigned int k=0;k<nq;++k)
- v_coarse(k*nd+d) = coarse.shape_value_component(face_c_dofs[i],k,d);
+ v_coarse(k*nd+d) = coarse.shape_value_component (face_c_dofs[i],k,d);
// solve the least squares
// problem.
const unsigned int nd = fe.n_components();
const unsigned int degree = fe.degree;
+ // prepare FEValues, quadrature etc on
+ // coarse cell
+ MappingCartesian<dim> mapping;
+ QGauss<dim> q_fine(degree+1);
+ const unsigned int nq = q_fine.size();
+
+ // create mass matrix on coarse cell.
+ FullMatrix<number> mass(n, n);
+ {
+ // set up a triangulation for coarse cell
+ Triangulation<dim,spacedim> tr;
+ GridGenerator::hyper_cube (tr, 0, 1);
+
+ FEValues<dim> coarse (mapping, fe, q_fine,
+ update_JxW_values | update_values);
+
+ typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
+ = tr.begin(0);
+ coarse.reinit (coarse_cell);
+
+ const std::vector<double> & JxW = coarse.get_JxW_values();
+ for (unsigned int i=0;i<n;++i)
+ for (unsigned int j=0;j<n;++j)
+ if (fe.is_primitive())
+ {
+ const double * coarse_i = &coarse.shape_value(i,0);
+ const double * coarse_j = &coarse.shape_value(j,0);
+ double mass_ij = 0;
+ for (unsigned int k=0;k<nq;++k)
+ mass_ij += JxW[k] * coarse_i[k] * coarse_j[k];
+ mass(i,j) = mass_ij;
+ }
+ else
+ {
+ double mass_ij = 0;
+ for (unsigned int d=0;d<nd;++d)
+ for (unsigned int k=0;k<nq;++k)
+ mass_ij += JxW[k] * coarse.shape_value_component(i,k,d)
+ * coarse.shape_value_component(j,k,d);
+ mass(i,j) = mass_ij;
+ }
+
+ // invert mass matrix
+ mass.gauss_jordan();
+ }
+
// loop over all possible refinement cases
for (unsigned int ref_case = RefinementCase<dim>::cut_x;
- ref_case<RefinementCase<dim>::isotropic_refinement+1; ++ref_case)
+ ref_case < RefinementCase<dim>::isotropic_refinement+1; ++ref_case)
{
const unsigned int
nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
Assert(matrices[ref_case-1][i].m() == n,
ExcDimensionMismatch(matrices[ref_case-1][i].m(),n));
}
-
+
+ // create a respective refinement on the
+ // triangulation
Triangulation<dim,spacedim> tr;
GridGenerator::hyper_cube (tr, 0, 1);
tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
tr.execute_coarsening_and_refinement();
-
- MappingCartesian<dim> mapping;
- QGauss<dim> q_fine(degree+1);
- const unsigned int nq = q_fine.size();
-
- FEValues<dim> coarse (mapping, fe, q_fine,
- update_quadrature_points |
- update_JxW_values |
- update_values);
+
FEValues<dim> fine (mapping, fe, q_fine,
- update_quadrature_points |
- update_JxW_values |
+ update_quadrature_points | update_JxW_values |
update_values);
-
+
typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
= tr.begin(0);
- // Compute the coarse level mass
- // matrix
- coarse.reinit(coarse_cell);
- FullMatrix<number> A(n, n);
- for (unsigned int k=0;k<nq;++k)
- for (unsigned int i=0;i<n;++i)
- for (unsigned int j=0;j<n;++j)
- if (fe.is_primitive())
- A(i,j) += coarse.JxW(k)
- * coarse.shape_value(i,k)
- * coarse.shape_value(j,k);
- else
- for (unsigned int d=0;d<nd;++d)
- A(i,j) = coarse.JxW(k)
- * coarse.shape_value_component(i,k,d)
- * coarse.shape_value_component(j,k,d);
-
- Householder<double> H(A);
Vector<number> v_coarse(n);
Vector<number> v_fine(n);
// Build RHS
+ const std::vector<double> & JxW = fine.get_JxW_values();
+
// Outer loop over all fine
// grid shape functions phi_j
for (unsigned int j=0;j<fe.dofs_per_cell;++j)
{
- v_fine = 0.;
- // Loop over all quadrature points
- for (unsigned int k=0;k<fine.n_quadrature_points;++k)
+ for (unsigned int i=0; i<fe.dofs_per_cell;++i)
{
- // integrate the scalar
- // product
- // (phi_i,phi_j) for
- // all coarse shape
- // functions to get the
- // right hand side
- for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+ if (fe.is_primitive())
{
- if (fe.is_primitive())
- v_fine(i) += fine.JxW(k)
- * coarse.shape_value(i,k)
- * fine.shape_value(j,k);
- else
- for (unsigned int d=0;d<nd;++d)
- v_fine(i) += fine.JxW(k)
- * coarse.shape_value_component(i,k,d)
- * fine.shape_value_component(j,k,d);
+ const double * coarse_i = &coarse.shape_value(i,0);
+ const double * fine_j = &fine.shape_value(j,0);
+
+ double update = 0;
+ for (unsigned int k=0; k<nq; ++k)
+ update += JxW[k] * coarse_i[k] * fine_j[k];
+ v_fine(i) = update;
+ }
+ else
+ {
+ double update = 0;
+ for (unsigned int d=0; d<nd; ++d)
+ for (unsigned int k=0; k<nq; ++k)
+ update += JxW[k] * coarse.shape_value_component(i,k,d)
+ * fine.shape_value_component(j,k,d);
+ v_fine(i) = update;
}
}
+
// RHS ready. Solve system
// and enter row into
// matrix
- H.least_squares(v_coarse, v_fine);
+ mass.vmult (v_coarse, v_fine);
for (unsigned int i=0;i<fe.dofs_per_cell;++i)
this_matrix(i,j) = v_coarse(i);
}
#include <base/thread_management.h>
#include <base/work_stream.h>
#include <base/multithread_info.h>
-#include <grid/tria_iterator.h>
-#include <grid/geometry_info.h>
+#include <base/geometry_info.h>
+#include <base/quadrature.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
-#include <grid/tria_iterator.h>
-#include <base/geometry_info.h>
-#include <base/quadrature.h>
#include <fe/fe.h>
#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+#include <grid/tria_iterator.h>
+#include <grid/geometry_info.h>
#include <hp/fe_values.h>
#include <hp/mapping_collection.h>
#include <numerics/matrices.h>
#include <lac/block_vector.h>
#include <lac/sparse_matrix.h>
#include <lac/block_sparse_matrix.h>
-#include <fe/mapping_q1.h>
#include <algorithm>
#include <set>
{
Scratch (const FiniteElement<dim,spacedim> &fe,
const UpdateFlags update_flags,
- const Function<spacedim> *coefficient,
- const Function<spacedim> *rhs_function,
+ const Function<spacedim> *coefficient,
+ const Function<spacedim> *rhs_function,
const Quadrature<dim> &quadrature,
- const Mapping<dim,spacedim> &mapping)
+ const Mapping<dim,spacedim> &mapping)
:
fe_collection (fe),
quadrature_collection (quadrature),
Scratch (const ::dealii::hp::FECollection<dim,spacedim> &fe,
const UpdateFlags update_flags,
- const Function<spacedim> *coefficient,
- const Function<spacedim> *rhs_function,
- const ::dealii::hp::QCollection<dim> &quadrature,
- const ::dealii::hp::MappingCollection<dim,spacedim> &mapping)
+ const Function<spacedim> *coefficient,
+ const Function<spacedim> *rhs_function,
+ const ::dealii::hp::QCollection<dim> &quadrature,
+ const ::dealii::hp::MappingCollection<dim,spacedim> &mapping)
:
fe_collection (fe),
quadrature_collection (quadrature),
update_flags (data.update_flags)
{}
- const ::dealii::hp::FECollection<dim,spacedim> fe_collection;
- const ::dealii::hp::QCollection<dim> quadrature_collection;
+ const ::dealii::hp::FECollection<dim,spacedim> fe_collection;
+ const ::dealii::hp::QCollection<dim> quadrature_collection;
const ::dealii::hp::MappingCollection<dim,spacedim> mapping_collection;
::dealii::hp::FEValues<dim,spacedim> x_fe_values;
std::vector<double> rhs_values;
std::vector<dealii::Vector<double> > rhs_vector_values;
+ std::vector<double> old_JxW;
+
const Function<spacedim> *coefficient;
const Function<spacedim> *rhs_function;
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
n_q_points = fe_values.n_quadrature_points;
- const FiniteElement<dim,spacedim> &fe = fe_values.get_fe();
+ const FiniteElement<dim,spacedim> &fe = fe_values.get_fe();
const unsigned int n_components = fe.n_components();
Assert(data.rhs_function == 0 ||
::dealii::MatrixCreator::ExcComponentMismatch());
copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
- copy_data.cell_matrix = 0;
-
copy_data.cell_rhs.reinit (dofs_per_cell);
- copy_data.cell_rhs = 0;
copy_data.dof_indices.resize (dofs_per_cell);
cell->get_dof_indices (copy_data.dof_indices);
- if (data.rhs_function != 0)
+ const bool use_rhs_function = data.rhs_function != 0;
+ if (use_rhs_function)
{
if (data.rhs_function->n_components==1)
{
}
}
- if (data.coefficient != 0)
+ const bool use_coefficient = data.coefficient != 0;
+ if (use_coefficient)
{
if (data.coefficient->n_components==1)
{
}
- if (data.coefficient != 0)
- {
- if (data.coefficient->n_components==1)
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ double add_data;
+ const std::vector<double> & JxW = fe_values.get_JxW_values();
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (fe.is_primitive ())
+ {
+ const unsigned int component_i =
+ fe.system_to_component_index(i).first;
+ const double * phi_i = &fe_values.shape_value(i,0);
+ add_data = 0;
+
+ // use symmetry in the mass matrix here:
+ // just need to calculate the diagonal
+ // and half of the elements above the
+ // diagonal
+ for (unsigned int j=i; j<dofs_per_cell; ++j)
+ if ((n_components==1) ||
+ (fe.system_to_component_index(j).first ==
+ component_i))
{
- const unsigned int component_i =
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first ==
- component_i))
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j)
- += (fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point) *
- fe_values.JxW(point) *
- data.coefficient_values[point]);
-
- if (data.rhs_function != 0)
+ const double * phi_j = &fe_values.shape_value(j,0);
+ add_data = 0;
+ if (use_coefficient)
{
- if (data.rhs_function->n_components==1)
+ if (data.coefficient->n_components==1)
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_values[point] * fe_values.JxW(point);
+ add_data += (phi_i[point] * phi_j[point] * JxW[point] *
+ data.coefficient_values[point]);
else
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_vector_values[point](component_i) *
- fe_values.JxW(point);
+ add_data += (phi_i[point] * phi_j[point] * JxW[point] *
+ data.coefficient_vector_values[point](component_i));
}
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ add_data += phi_i[point] * phi_j[point] * JxW[point];
+
+ // this is even ok for i==j, since then
+ // we just write the same value twice.
+ copy_data.cell_matrix(i,j) = add_data;
+ copy_data.cell_matrix(j,i) = add_data;
}
- }
- else
- {
- if (fe.is_primitive ())
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+
+ if (use_rhs_function)
+ {
+ add_data = 0;
+ if (data.rhs_function->n_components==1)
+ for (unsigned int point=0; point<n_q_points; ++point)
+ add_data += phi_i[point] * JxW[point] *
+ data.rhs_values[point];
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ add_data += phi_i[point] * JxW[point] *
+ data.rhs_vector_values[point](component_i);
+ copy_data.cell_rhs(i) = add_data;
+ }
+ }
+ else
+ {
+ // non-primitive vector-valued FE, using
+ // symmetry again
+ for (unsigned int j=i; j<dofs_per_cell; ++j)
+ {
+ add_data = 0;
+ for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+ if (fe.get_nonzero_components(i)[comp_i] &&
+ fe.get_nonzero_components(j)[comp_i])
{
- const unsigned int component_i =
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first ==
- component_i))
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j) +=
- (fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point) *
- fe_values.JxW(point) *
- data.coefficient_vector_values[point](component_i));
-
- if (data.rhs_function != 0)
+ if (use_coefficient)
{
- if (data.rhs_function->n_components==1)
+ if (data.coefficient->n_components==1)
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_values[point] * fe_values.JxW(point);
+ add_data += (fe_values.shape_value_component(i,point,comp_i) *
+ fe_values.shape_value_component(j,point,comp_i) *
+ JxW[point] *
+ data.coefficient_values[point]);
else
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_vector_values[point](component_i) *
- fe_values.JxW(point);
+ add_data += (fe_values.shape_value_component(i,point,comp_i) *
+ fe_values.shape_value_component(j,point,comp_i) *
+ JxW[point] *
+ data.coefficient_vector_values[point](comp_i));
}
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ add_data += fe_values.shape_value_component(i,point,comp_i) *
+ fe_values.shape_value_component(j,point,comp_i) * JxW[point];
}
- }
- else
- // non-primitive vector-valued FE
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i])
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
- if (fe.get_nonzero_components(j)[comp_j])
- if (comp_i == comp_j)
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j) +=
- (fe_values.shape_value_component(i,point,comp_i) *
- fe_values.shape_value_component(j,point,comp_j) *
- fe_values.JxW(point) *
- data.coefficient_vector_values[point](comp_i));
-
- if (data.rhs_function != 0)
- {
- if (data.rhs_function->n_components==1)
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) +=
- fe_values.shape_value_component(i,point,comp_i) *
- data.rhs_values[point] * fe_values.JxW(point);
- else
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) +=
- fe_values.shape_value_component(i,point,comp_i) *
- data.rhs_vector_values[point](comp_i) *
- fe_values.JxW(point);
- }
- }
- }
- }
- }
- else
- // no coefficient
- {
- if (fe.is_primitive ())
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const unsigned int component_i =
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first ==
- component_i))
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j) +=
- (fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point) *
- fe_values.JxW(point));
-
- if (data.rhs_function != 0)
+
+ copy_data.cell_matrix(i,j) = add_data;
+ copy_data.cell_matrix(j,i) = add_data;
+ }
+
+ if (use_rhs_function)
+ {
+ add_data = 0;
+ for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+ if (fe.get_nonzero_components(i)[comp_i])
{
if (data.rhs_function->n_components==1)
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_values[point] * fe_values.JxW(point);
+ add_data += fe_values.shape_value_component(i,point,comp_i) *
+ JxW[point] * data.rhs_values[point];
else
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_vector_values[point](component_i) *
- fe_values.JxW(point);
- }
- }
- }
- else
- // non-primitive FE, no coefficient
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i])
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
- if (fe.get_nonzero_components(j)[comp_j])
- if (comp_i == comp_j)
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j) +=
- (fe_values.shape_value_component(i,point,comp_i) *
- fe_values.shape_value_component(j,point,comp_j) *
- fe_values.JxW(point));
-
- if (data.rhs_function != 0)
- {
- if (data.rhs_function->n_components==1)
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) +=
- fe_values.shape_value_component(i,point,comp_i) *
- data.rhs_values[point] * fe_values.JxW(point);
- else
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) +=
- fe_values.shape_value_component(i,point,comp_i) *
- data.rhs_vector_values[point](comp_i) *
- fe_values.JxW(point);
- }
+ add_data += fe_values.shape_value_component(i,point,comp_i) *
+ JxW[point] * data.rhs_vector_values[point](comp_i);
}
- }
- }
+ copy_data.cell_rhs(i) = add_data;
+ }
+ }
}
::dealii::MatrixCreator::ExcComponentMismatch());
copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
- copy_data.cell_matrix = 0;
-
copy_data.cell_rhs.reinit (dofs_per_cell);
- copy_data.cell_rhs = 0;
-
copy_data.dof_indices.resize (dofs_per_cell);
cell->get_dof_indices (copy_data.dof_indices);
- if (data.rhs_function != 0)
+ const bool use_rhs_function = data.rhs_function != 0;
+ if (use_rhs_function)
{
if (data.rhs_function->n_components==1)
{
}
}
- if (data.coefficient != 0)
+ const bool use_coefficient = data.coefficient != 0;
+ if (use_coefficient)
{
if (data.coefficient->n_components==1)
{
}
- if (data.coefficient != 0)
- {
- if (data.coefficient->n_components==1)
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ const std::vector<double> & JxW = fe_values.get_JxW_values();
+ double add_data;
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (fe.is_primitive ())
+ {
+ const unsigned int component_i =
+ fe.system_to_component_index(i).first;
+ const Tensor<1,spacedim> * grad_phi_i =
+ &fe_values.shape_grad(i,0);
+
+ // can use symmetry
+ for (unsigned int j=i; j<dofs_per_cell; ++j)
+ if ((n_components==1) ||
+ (fe.system_to_component_index(j).first ==
+ component_i))
{
- const unsigned int component_i =
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first ==
- component_i))
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j)
- += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) *
- fe_values.JxW(point) *
- data.coefficient_values[point]);
-
- if (data.rhs_function != 0)
+ const Tensor<1,spacedim> * grad_phi_j =
+ & fe_values.shape_grad(j,0);
+ add_data = 0;
+ if (use_coefficient)
{
- if (data.rhs_function->n_components==1)
+ if (data.coefficient->n_components==1)
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_values[point] * fe_values.JxW(point);
+ add_data += ((grad_phi_i[point]*grad_phi_j[point]) *
+ JxW[point] *
+ data.coefficient_values[point]);
else
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_vector_values[point](component_i) *
- fe_values.JxW(point);
+ add_data += ((grad_phi_i[point]*grad_phi_j[point]) *
+ JxW[point] *
+ data.coefficient_vector_values[point](component_i));
}
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ add_data += (grad_phi_i[point]*grad_phi_j[point]) *
+ JxW[point];
+
+ copy_data.cell_matrix(i,j) = add_data;
+ copy_data.cell_matrix(j,i) = add_data;
}
- }
- else
- {
- if (fe.is_primitive ())
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+
+ if (use_rhs_function)
+ {
+ const double * phi_i = &fe_values.shape_value(i,0);
+ add_data = 0;
+ if (data.rhs_function->n_components==1)
+ for (unsigned int point=0; point<n_q_points; ++point)
+ add_data += phi_i[point] * JxW[point] *
+ data.rhs_values[point];
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ add_data += phi_i[point] * JxW[point] *
+ data.rhs_vector_values[point](component_i);
+ copy_data.cell_rhs(i) = add_data;
+ }
+ }
+ else
+ {
+ // non-primitive vector-valued FE
+ for (unsigned int j=i; j<dofs_per_cell; ++j)
+ {
+ add_data = 0;
+ for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+ if (fe.get_nonzero_components(i)[comp_i] &&
+ fe.get_nonzero_components(j)[comp_i])
{
- const unsigned int component_i =
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first ==
- component_i))
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j) +=
- (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) *
- fe_values.JxW(point) *
- data.coefficient_vector_values[point](component_i));
-
- if (data.rhs_function != 0)
+ if (use_coefficient)
{
- if (data.rhs_function->n_components==1)
+ if (data.coefficient->n_components==1)
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_values[point] * fe_values.JxW(point);
+ add_data += ((fe_values.shape_grad_component(i,point,comp_i) *
+ fe_values.shape_grad_component(j,point,comp_i)) *
+ JxW[point] *
+ data.coefficient_values[point]);
else
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_vector_values[point](component_i) *
- fe_values.JxW(point);
+ add_data += ((fe_values.shape_grad_component(i,point,comp_i) *
+ fe_values.shape_grad_component(j,point,comp_i)) *
+ JxW[point] *
+ data.coefficient_vector_values[point](comp_i));
}
+ else
+ for (unsigned int point=0; point<n_q_points; ++point)
+ add_data += (fe_values.shape_grad_component(i,point,comp_i) *
+ fe_values.shape_grad_component(j,point,comp_i)) *
+ JxW[point];
}
- }
- else
- // non-primitive vector-valued FE
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i])
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
- if (fe.get_nonzero_components(j)[comp_j])
- if (comp_i == comp_j)
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j) +=
- (fe_values.shape_grad_component(i,point,comp_i) *
- fe_values.shape_grad_component(j,point,comp_j) *
- fe_values.JxW(point) *
- data.coefficient_vector_values[point](comp_i));
-
- if (data.rhs_function != 0)
- {
- if (data.rhs_function->n_components==1)
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) +=
- fe_values.shape_value_component(i,point,comp_i) *
- data.rhs_values[point] * fe_values.JxW(point);
- else
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) +=
- fe_values.shape_value_component(i,point,comp_i) *
- data.rhs_vector_values[point](comp_i) *
- fe_values.JxW(point);
- }
- }
- }
- }
- }
- else
- // no coefficient
- {
- if (fe.is_primitive ())
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const unsigned int component_i =
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first ==
- component_i))
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j) +=
- (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) *
- fe_values.JxW(point));
-
- if (data.rhs_function != 0)
+
+ copy_data.cell_matrix(i,j) = add_data;
+ copy_data.cell_matrix(j,i) = add_data;
+ }
+
+ if (use_rhs_function)
+ {
+ add_data = 0;
+ for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+ if (fe.get_nonzero_components(i)[comp_i])
{
if (data.rhs_function->n_components==1)
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_values[point] * fe_values.JxW(point);
+ add_data += fe_values.shape_value_component(i,point,comp_i) *
+ JxW[point] * data.rhs_values[point];
else
for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) += fe_values.shape_value(i, point) *
- data.rhs_vector_values[point](component_i) *
- fe_values.JxW(point);
- }
- }
- }
- else
- // non-primitive FE, no coefficient
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i])
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- for (unsigned int comp_j = 0; comp_j < n_components; ++comp_j)
- if (fe.get_nonzero_components(j)[comp_j])
- if (comp_i == comp_j)
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_matrix(i,j) +=
- (fe_values.shape_grad_component(i,point,comp_i) *
- fe_values.shape_grad_component(j,point,comp_j) *
- fe_values.JxW(point));
-
- if (data.rhs_function != 0)
- {
- if (data.rhs_function->n_components==1)
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) +=
- fe_values.shape_value_component(i,point,comp_i) *
- data.rhs_values[point] * fe_values.JxW(point);
- else
- for (unsigned int point=0; point<n_q_points; ++point)
- copy_data.cell_rhs(i) +=
- fe_values.shape_value_component(i,point,comp_i) *
- data.rhs_vector_values[point](comp_i) *
- fe_values.JxW(point);
- }
+ add_data += fe_values.shape_value_component(i,point,comp_i) *
+ JxW[point] * data.rhs_vector_values[point](comp_i);
}
- }
- }
+ copy_data.cell_rhs(i) = add_data;
+ }
+ }
}
-
-
+
+
template <typename MatrixType,
typename VectorType>
void copy_local_to_global (const AssemblerData::CopyData &data,
(data.cell_rhs.size() == dofs_per_cell),
ExcInternalError());
- matrix->add (data.dof_indices, data.cell_matrix);
-
+ matrix->add(data.dof_indices, data.cell_matrix);
if (right_hand_side != 0)
for (unsigned int i=0; i<dofs_per_cell; ++i)
(*right_hand_side)(data.dof_indices[i]) += data.cell_rhs(i);