From: goll Date: Thu, 16 Feb 2012 15:16:54 +0000 (+0000) Subject: Added the possibility to compute laplacians of a FEFieldFunction. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=797c884ac834f697da71af422de0f2cc0e9006a1;p=dealii-svn.git Added the possibility to compute laplacians of a FEFieldFunction. git-svn-id: https://svn.dealii.org/trunk@25100 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/numerics/fe_field_function.h b/deal.II/include/deal.II/numerics/fe_field_function.h index 151c5255e9..fcb22bdfee 100644 --- a/deal.II/include/deal.II/numerics/fe_field_function.h +++ b/deal.II/include/deal.II/numerics/fe_field_function.h @@ -297,6 +297,41 @@ namespace Functions gradient_list (const std::vector< Point< dim > > &p, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const; + + + /** + * Compute the Laplacian of a + * given component at point p. + */ + virtual double + laplacian (const Point &p, + const unsigned int component = 0) const; + + /** + * Compute the Laplacian of all + * components at point p and + * store them in values. + */ + virtual void + vector_laplacian (const Point &p, + Vector &values) const; + + /** + * Compute the Laplacian of one + * component at a set of points. + */ + virtual void + laplacian_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; + + /** + * Compute the Laplacians of all + * components at a set of points. + */ + virtual void + vector_laplacian_list (const std::vector > &points, + std::vector > &values) const; /** * Create quadrature diff --git a/deal.II/include/deal.II/numerics/fe_field_function.templates.h b/deal.II/include/deal.II/numerics/fe_field_function.templates.h index 70dce9ed44..a5687b4e28 100644 --- a/deal.II/include/deal.II/numerics/fe_field_function.templates.h +++ b/deal.II/include/deal.II/numerics/fe_field_function.templates.h @@ -139,6 +139,52 @@ namespace Functions return grads[comp]; } + + + template + void FEFieldFunction::vector_laplacian + (const Point &p,Vector &values) const + { + Assert (values.size() == n_components, + ExcDimensionMismatch(values.size(), n_components)); + typename DH::active_cell_iterator cell = cell_hint.get(); + if (cell == dh->end()) + cell = dh->begin_active(); + Point qp = mapping.transform_real_to_unit_cell(cell, p); + + // Check if we already have all we need + if (!GeometryInfo::is_inside_unit_cell(qp)) + { + const std::pair > my_pair + = GridTools::find_active_cell_around_point (mapping, *dh, p); + cell = my_pair.first; + qp = my_pair.second; + } + + cell_hint.get() = cell; + + // Now we can find out about the point + Quadrature quad(qp); + FEValues fe_v(mapping, dh->get_fe(), quad, + update_hessians); + fe_v.reinit(cell); + std::vector< Vector > vvalues (1, values); + fe_v.get_function_laplacians(data_vector, vvalues); + values = vvalues[0]; + } + + + + template + double FEFieldFunction::laplacian + (const Point &p, const unsigned int comp) const + { + Vector lap(n_components); + vector_laplacian(p, lap); + return lap[comp]; + } + + // Now the list versions // ============================== @@ -249,6 +295,58 @@ namespace Functions values[q] = vvalues[q][component]; } + + template + void + FEFieldFunction:: + vector_laplacian_list (const std::vector > & points, + std::vector< Vector > &values) const + { + Assert(points.size() == values.size(), + ExcDimensionMismatch(points.size(), values.size())); + + std::vector cells; + std::vector > > qpoints; + std::vector > maps; + + unsigned int ncells = compute_point_locations(points, cells, qpoints, maps); + + // Now gather all the informations we need + for (unsigned int i=0; i ww(nq, 1./((double) nq)); + Quadrature quad(qpoints[i], ww); + + // Get a function value object + FEValues fe_v(mapping, dh->get_fe(), quad, + update_hessians); + fe_v.reinit(cells[i]); + std::vector< Vector > vvalues (nq, Vector(n_components)); + fe_v.get_function_laplacians(data_vector, vvalues); + for (unsigned int q=0; q + void + FEFieldFunction:: + laplacian_list (const std::vector > &points, + std::vector< double > &values, + const unsigned int component) const + { + Assert(points.size() == values.size(), + ExcDimensionMismatch(points.size(), values.size())); + std::vector< Vector > vvalues(points.size(), Vector(n_components)); + vector_laplacian_list(points, vvalues); + for (unsigned int q=0; q