From: hartmann Date: Thu, 1 Jul 1999 19:56:15 +0000 (+0000) Subject: implementation of #interpolate# and #integrate_difference# for vector functions X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=80b64d4163a9a295adfd44cc6eabb3bd13e7713d;p=dealii-svn.git implementation of #interpolate# and #integrate_difference# for vector functions git-svn-id: https://svn.dealii.org/trunk@1525 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 87154e5dd1..9ed65f0b1b 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -73,6 +73,202 @@ void VectorTools::interpolate (const DoFHandler &dof, }; + +template +void VectorTools::interpolate (const DoFHandler &dof, + const VectorFunction&vectorfunction, + Vector &vec) +{ + const FiniteElement &fe = dof.get_fe(); + + DoFHandler::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + + // For FESystems many of the + // unit_support_points will + // appear multiply, as a point + // may be unit_support_point + // for several of the components + // of the system. + // The following is rather + // complicated as it is + // avoided to evaluate + // the vectorfunction multiply at + // the same point on a cell. + vector > unit_support_points (fe.total_dofs); + fe.get_unit_support_points(unit_support_points); + + // The following works well + // if #dofs_per_cell<=1# as then + // the multiple support_points + // are placed one after another. + + // find the support points + // on a cell that + // are multiply mentioned in + // #unit_support_points#. + // Mark the first representative + // of each multiply mentioned + // support point by setting + // #true# in the boolean vector + // #is_representative_point#. +// vector is_representative_point(fe.total_dofs, false); +// is_representative_point[0]=true; +// unsigned int n_rep_points=1; +// for (unsigned int last_rep_point=0, i=1; i dofs_on_cell (fe.total_dofs); +// vector > support_points (fe.total_dofs); + +// vector > rep_points (n_rep_points); +// vector > function_values_at_rep_points ( +// n_rep_points, Vector(fe.n_components)); + +// for (; cell!=endc; ++cell) +// { +// // for each cell: +// // get location of finite element +// // off-points (support_points) +// fe.get_support_points (cell, support_points); + +// // pick out the representative +// // support points +// unsigned int j=0; +// for (unsigned int i=0; iget_dof_indices (dofs_on_cell); + +// // distribute function values to the +// // whole vector +// int last_rep_point = -1; +// // it holds `is_representative_point[0]=true' +// // therefore the first #last_rep_point# is 0 +// // and we need to start with +// // `last_rep_point = -1' +// for (unsigned int i=0; i1#, + // i.e. it is usable also for systems + // including + // FEQ3, FEQ4, FEDG_Qx. + + // Find the support points + // on a cell that + // are multiply mentioned in + // #unit_support_points#. + // Mark the first representative + // of each multiply mentioned + // support point by appending its + // dof index to #dofs_of_rep_points#. + // Each multiple point gets to know + // the dof index of its representative + // point by the #dof_to_rep_dof_table#. + + // the following vector collects all dofs i, + // 0<=i dofs_of_rep_points; + // the following table converts a dof i + // to the dof of the representative + // point. + vector dof_to_rep_dof_table; + unsigned int n_rep_points=0; + for (unsigned int i=0; i0; --j) + if (unit_support_points[i] + == unit_support_points[dofs_of_rep_points[j-1]]) + { + dof_to_rep_dof_table.push_back(j-1); + representative=false; + break; + } + + if (representative) + { + dofs_of_rep_points.push_back(i); + dof_to_rep_dof_table.push_back(i); + ++n_rep_points; + } + } + Assert(dofs_of_rep_points.size()==n_rep_points, ExcInternalError()); + Assert(dof_to_rep_dof_table.size()==fe.total_dofs, ExcInternalError()); + + cout << "n_rep_points=" << n_rep_points << endl; + + vector dofs_on_cell (fe.total_dofs); + vector > support_points (fe.total_dofs); + + vector > rep_points (n_rep_points); + vector > function_values_at_rep_points ( + n_rep_points, Vector(fe.n_components)); + + for (; cell!=endc; ++cell) + { + // for each cell: + // get location of finite element + // off-points (support_points) + fe.get_support_points (cell, support_points); + + // pick out the representative + // support points + for (unsigned int j=0; jget_dof_indices (dofs_on_cell); + + // distribute the function values to + // the global vector + for (unsigned int i=0; i void VectorTools::interpolate(const DoFHandler &high_dof, const DoFHandler &low_dof, @@ -684,6 +880,8 @@ VectorTools::integrate_difference (const DoFHandler &dof, const Quadrature &q, const NormType &norm) { + Assert(norm != mean , ExcNotUseful()); + const FiniteElement &fe = dof.get_fe(); difference.reinit (dof.get_tria().n_active_cells()); @@ -703,102 +901,89 @@ VectorTools::integrate_difference (const DoFHandler &dof, // initialize for this cell fe_values.reinit (cell); - switch (norm) + switch (norm) { case mean: + break; 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 > psi (n_q_points); - // in praxi: first compute - // exact fe_function vector - exact_solution.value_list (fe_values.get_quadrature_points(), - psi); + // first compute the exact solution + // (vectors) at the quadrature points + exact_solution.value_list (fe_values.get_quadrature_points(), psi); // then subtract finite element // fe_function if (true) { - vector< Vector > function_values (n_q_points, - Vector(fe.n_components)); + vector< Vector > function_values ( + n_q_points, Vector(fe.n_components)); + fe_values.get_function_values (fe_function, function_values); -/* transform (psi.begin(), psi.end(), - function_values.begin(), - psi.begin(), - minus()); -*/ }; + for (unsigned int q=0; q psi_scalar (n_q_points); + 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)); + for (unsigned int q=0; q::integrate_difference (const DoFHandler &dof, // H1_norm starts at the previous // case statement, but continues // here! - - // for H1_norm: re-square L2_norm. - diff = sqr(diff); + // Until now, #diff# includes the + // square of the L2_norm. // same procedure as above, but now - // psi is a vector of gradients + // psi is a vector of Jacobians + // i.e. psi is a vector of vectors of + // gradients. const unsigned int n_q_points = q.n_quadrature_points; - vector > psi (n_q_points); - + vector > > psi ( + n_q_points, vector >(fe.n_components, Tensor<1,dim>())); + // in praxi: first compute // exact fe_function vector -/* exact_solution.gradient_list (fe_values.get_quadrature_points(), - psi); -*/ + exact_solution.gradient_list (fe_values.get_quadrature_points(), psi); + // then subtract finite element - // fe_function + // function_grads if (true) { - vector > function_grads (n_q_points, Tensor<1,dim>()); + vector > > function_grads ( + n_q_points, vector >(fe.n_components, Tensor<1,dim>())); fe_values.get_function_grads (fe_function, function_grads); -/* transform (psi.begin(), psi.end(), - function_grads.begin(), - psi.begin(), - minus >()); -*/ }; + for (unsigned int q=0; q psi_square (psi.size(), 0.0); - for (unsigned int i=0; i