From be47c5f02e608bc74f4ed53b90e4b486e956f2ec Mon Sep 17 00:00:00 2001 From: heltai Date: Wed, 29 Aug 2007 20:16:16 +0000 Subject: [PATCH] Added FEFieldFunction and ParsedFunction classes. git-svn-id: https://svn.dealii.org/trunk@15091 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/parsed_function.h | 197 +++++++++++ deal.II/base/source/parsed_function.cc | 122 +++++++ .../include/numerics/fe_field_function.h | 250 +++++++++++++ .../numerics/fe_field_function.templates.h | 327 ++++++++++++++++++ .../source/numerics/fe_field_function.cc | 76 ++++ deal.II/doc/authors.html | 4 +- deal.II/doc/news/changes.html | 13 + 7 files changed, 987 insertions(+), 2 deletions(-) create mode 100644 deal.II/base/include/base/parsed_function.h create mode 100644 deal.II/base/source/parsed_function.cc create mode 100644 deal.II/deal.II/include/numerics/fe_field_function.h create mode 100644 deal.II/deal.II/include/numerics/fe_field_function.templates.h create mode 100644 deal.II/deal.II/source/numerics/fe_field_function.cc diff --git a/deal.II/base/include/base/parsed_function.h b/deal.II/base/include/base/parsed_function.h new file mode 100644 index 0000000000..69da46e8dc --- /dev/null +++ b/deal.II/base/include/base/parsed_function.h @@ -0,0 +1,197 @@ +//--------------------------------------------------------------------------- +// $Id: function_parser.h 14594 2007-03-22 20:17:41Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2007 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#ifndef __deal2__parsed_function_h +#define __deal2__parsed_function_h + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions { + /** Friendly interface to the FunctionParser class. This class is + meant as a wrapper for the FunctionParser class. It provides two + methods to declare and parse a ParameterHandler object and creates + the Function object declared in the parameter file. This class is + derived from the AutoDerivativeFunction class, so you don't need + to specify derivatives. An example of usage of this class is as follows: + + \code + // A parameter handler + ParameterHandler prm; + + // Declare a section for the function we need + prm.enter_subsection("My vector function"); + ParsedFunction::declare_parameters(prm, dim); + prm.leave_subsection(); + + // Create a ParsedFunction + ParsedFunction my_vector_function(dim); + + // Parse an input file. + prm.read_input(some_input_file); + + // Initialize the ParsedFunction object with the given file + prm.enter_subsection("My vector function"); + my_vector_function.parse_parameters(prm); + prm.leave_subsection(); + + \endcode + + And here is an example of how the input parameter could look like + (see the documentation of the FunctionParser class for a detailed + description of the syntax of the function definition): + + \code + + # A test two dimensional vector function, depending on time + subsection My vector function + set Function constants = kappa=.1, lambda=2. + set Function expression = if(y>.5, kappa*x*(1-x),0); t^2*cos(lambda*pi*x) + set Variable names = x,y,t + end + + \endcode + + \ingroup functions + \author Luca Heltai, 2006 + */ + template + class ParsedFunction : public AutoDerivativeFunction + { + public: + /** Construct a vector function. The vector function which is + generated has @p n_components components (defaults to 1). The parameter + @p h is used to initialize the AutoDerivativeFunction class from + which this class is derived. */ + ParsedFunction (const unsigned int n_components = 1, const double h=1e-8); + + /** Declare parameters needed by this class. The additional + parameter @p n_components is used to generate the right code according + to the number of components of the function that will parse this + ParameterHandler. If the number of components which is parsed + does not match the number of components of this object, an + assertion is thrown and the program is aborted. + + The default behavior for this class is to declare the following + entries: + + \code + + set Function constants = + set Function expression = 0 + set Variable names = x,y,t + + \endcode + + */ + static void declare_parameters(ParameterHandler &prm, + const unsigned int n_components = 1); + + /** Parse parameters needed by this class. If the number of + components which is parsed does not match the number of + components of this object, an assertion is thrown and the + program is aborted. + + In order for the class to function properly, we follow the same + convenctions declared in the FunctionParser class (look there + for a detailed description of the syntax for function + declarations). + + The three variables that can be parsed from a parameter file are + the following: + + \code + + set Function constants = + set Function expression = + set Variable names = + + \endcode + + Function constants is a collection of pairs in the form + name=value, separated by commas, for example: + + \code + + set Function constants = lambda=1. , alpha=2., gamma=3. + + \endcode + + These constants can be used in the declaration of the function + expression, which follows the convention of the FunctionParser + class. In order to specify vector functions, semicolons have to + be used to separate the different components, e.g.: + + \code + + set Function expression = cos(pi*x) ; cos(pi*y) + + \endcode + + The variable names entry can be used to customize the name of + the variables used in the Function. It defaults to + + \code + + set Variable names = x,t + + \endcode + + for one dimensional problems, + + \code + + set Variable names = x,y,t + + \endcode + + for two dimensional problems and + + \code + + set Variable names = x,y,z,t + + \endcode + + for three dimensional problems. + + The time variable can be set according to specifications in the + FunctionTime class. + + */ + void parse_parameters(ParameterHandler &prm); + + /** Get one value at the given point. */ + virtual void vector_value (const Point &p, + Vector &values) const; + + /** Return the value of the function at the given point. Unless + there is only one component (i.e. the function is scalar), you + should state the component you want to have evaluated; it + defaults to zero, i.e. the first component. */ + virtual double value (const Point< dim > & p, + const unsigned int component = 0) const; + + /** We need to overwrite this to set the time also in the accessor + FunctionParser. */ + virtual void set_time(const double newtime); + private: + FunctionParser function_object; + }; +} +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/base/source/parsed_function.cc b/deal.II/base/source/parsed_function.cc new file mode 100644 index 0000000000..995eb64bb6 --- /dev/null +++ b/deal.II/base/source/parsed_function.cc @@ -0,0 +1,122 @@ +//--------------------------------------------------------------------------- +// $Id: function_parser.h 14594 2007-03-22 20:17:41Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2007 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions { + template + ParsedFunction::ParsedFunction (const unsigned int n_components, const double h) : + AutoDerivativeFunction(h, n_components), + function_object(n_components) + { + } + + template + void ParsedFunction::declare_parameters(ParameterHandler &prm, unsigned int n_components) + { + Assert(n_components > 0, ExcZero()); + + std::string vnames; + switch (dim) { + case 1: + vnames = "x,t"; + break; + case 2: + vnames = "x,y,t"; + break; + case 3: + vnames = "x,y,z,t"; + break; + default: + AssertThrow(false, ExcNotImplemented()); + break; + } + prm.declare_entry("Variable names", vnames, Patterns::Anything(), + "The name of the variables as they will be used in the function, separated by ','."); + // The expression of the function + std::string expr = "0"; + for(unsigned int i=1; i + void ParsedFunction::parse_parameters(ParameterHandler &prm) + { + std::string vnames = prm.get("Variable names"); + std::string expression = prm.get("Function expression"); + std::string constants_list = prm.get("Function constants"); + + std::vector const_list = + Utilities::split_string_list(constants_list, ','); + std::map constants; + for(unsigned int i = 0; i < const_list.size(); ++i) { + std::vector this_c = + Utilities::split_string_list(const_list[i], '='); + AssertThrow(this_c.size() == 2, ExcMessage("Invalid format")); + double tmp; + AssertThrow( sscanf(this_c[1].c_str(), "%lf", &tmp), ExcMessage("Double number?")); + constants[this_c[0]] = tmp; + } + + constants["pi"] = deal_II_numbers::PI; + constants["Pi"] = deal_II_numbers::PI; + + unsigned int nn = (Utilities::split_string_list(vnames)).size(); + switch (nn) { + case dim: + // Time independent function + function_object.initialize(vnames, expression, constants); + break; + case dim+1: + // Time dependent function + function_object.initialize(vnames, expression, constants, true); + break; + default: + AssertThrow(false, ExcMessage("Not the correct size. Check your code.")); + } + } + + template + void ParsedFunction::vector_value (const Point &p, + Vector &values) const + { + function_object.vector_value(p, values); + } + + template + double ParsedFunction::value (const Point &p, + unsigned int comp) const + { + return function_object.value(p, comp); + } + + template + void ParsedFunction::set_time (const double newtime) + { + function_object.set_time(newtime); + AutoDerivativeFunction::set_time(newtime); + } + + // Explicit instantiations + template class ParsedFunction<1>; + template class ParsedFunction<2>; + template class ParsedFunction<3>; +} +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/include/numerics/fe_field_function.h b/deal.II/deal.II/include/numerics/fe_field_function.h new file mode 100644 index 0000000000..3137e18e31 --- /dev/null +++ b/deal.II/deal.II/include/numerics/fe_field_function.h @@ -0,0 +1,250 @@ +//--------------------------------------------------------------------------- +// $Id: function_parser.h 14594 2007-03-22 20:17:41Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2007 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__fe_function_h +#define __deal2__fe_function_h + +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions { + + /** This is an interpolation function for the given dof handler and + the given solution vector. The points at which this function can + be evaluated MUST be inside the domain of the dof handler, but + except from this, no other requirement is given. This function is + rather slow, as it needs to construct a quadrature object for the + point (or set of points) where you want to evaluate your finite + element function. In order to do so, it needs to find out where + the points lie. + + If you know in advance in which cell your points lye, you can + accelerate things a bit, by calling set_active_cell before + asking for values or gradients of the function. If you don't do + this, and your points don't lie in the cell that is currently + stored, the function GridTools::find_cell_around_point is called + to find out where the point is. You can specify an optional + mapping to use when looking for points in the grid. If you don't + do so, this function uses a Q1 mapping. + + Once the FEFieldFunction knows where the points lie, it creates a + quadrature formula for those points, and calls + FEValues::get_function_values or FEValues::get_function_grads with + the given quadrature points. + + If you only need the quadrature points but not the values of the + finite element function (you might want this for the adjoint + interpolation), you can also use the function @p + compute_point_locations alone. + + An example of how to use this function is the following: + + \code + + // Generate two triangulations + Triangulation tria_1; + Triangulation tria_2; + + // Read the triangulations from files, or build them up, or get + // them from some place... Assume that tria_2 is *entirely* + // included in tria_1 + ... + + // Associate a dofhandler and a solution to the first + // triangulation + DoFHandler dh1(tria_1); + Vector solution_1; + + // Do the same with the second + DoFHandler dh2; + Vector solution_2; + + // Setup the system, assemble matrices, solve problems and get the + // nobel prize on the first domain... + ... + + // Now project it to the second domain + FEFieldFunction fe_function_1 (dh_1, solution_1); + VectorTools::project(dh_2, constraints_2, quad, fe_function_1, solution_2); + + // Or interpolate it... + Vector solution_3; + VectorTools::interpolate(dh_2, fe_function_1, solution_3); + + \endcode + + The snippet of code above will work assuming that the second + triangulation is entirely included in the first one. + + FEFieldFunction is designed to be an easy way to get the results of + your computations across different, possibly non matching, + grids. No knowledge of the location of the points is assumed in + this class, which makes it rely entirely on the + GridTools::find_active_cell_around_point utility for its + job. However the class can be fed an "educated guess" of where the + points that will be computed actually are by using the + FEFieldFunction::set_active_cell method, so if you have a smart way to + tell where your points are, you will save a lot of computational + time by letting this class know. + + An optimization based on a caching mechanism was used by the + author of this class for the implementation of a Finite Element + Immersed Boundary Method. + + \addtogroup functions + + \author Luca Heltai, 2006 + + \todo Add hp functionality + */ + template , + typename VECTOR=Vector > + class FEFieldFunction : public Function + { + public: + /** Construct a vector function. A smart pointers is stored to the + dof handler, so you have to make sure that it make sense for + the entire lifetime of this object. The number of components + of this functions is equal to the number of components of the + finite element object. If a mapping is specified, that is what + is used to find out where the points lay. Otherwise the + standard Q1 mapping is used. */ + FEFieldFunction (const DH &dh, const VECTOR &data_vector, + const Mapping &mapping = StaticMappingQ1::mapping); + + /** Set the current cell. If you know in advance where your points + lie, you can tell this object by calling this function. This + will speed things up a little. */ + inline void set_active_cell(typename DH::active_cell_iterator &newcell); + + /** Get ONE vector value at the given point. It is inefficient to + use single points. If you need more than one at a time, use the + vector_value_list function. For efficiency reasons, it is better + if all the points lie on the same cell. This is not mandatory, + however it does speed things up. */ + virtual void vector_value (const Point &p, + Vector &values) const; + + /** Return the value of the function at the given point. Unless + there is only one component (i.e. the function is scalar), you + should state the component you want to have evaluated; it + defaults to zero, i.e. the first component. It is inefficient + to use single points. If you need more than one at a time, use + the vector_value_list function. For efficiency reasons, it is + better if all the points lie on the same cell. This is not + mandatory, however it does speed things up. */ + virtual double value (const Point< dim > & p, + const unsigned int component = 0) const; + + /** Set @p values to the point values of the specified component of + the function at the @p points. It is assumed that @p values + already has the right size, i.e. the same size as the points + array. This is rather efficient if all the points lie on the + same cell. If this is not the case, things may slow down a bit. + */ + virtual void value_list (const std::vector > & points, + std::vector< double > &values, + const unsigned int component = 0) const; + + + /** Set @p values to the point values of the function at the @p + points. It is assumed that @p values already has the right size, + i.e. the same size as the points array. This is rather efficient + if all the points lie on the same cell. If this is not the case, + things may slow down a bit. + */ + virtual void vector_value_list (const std::vector > & points, + std::vector< Vector > &values) const; + + /** Return the gradient of all components of the function at the + given point. It is inefficient to use single points. If you + need more than one at a time, use the vector_value_list + function. For efficiency reasons, it is better if all the points + lie on the same cell. This is not mandatory, however it does + speed things up. */ + virtual void + vector_gradient (const Point< dim > &p, + std::vector< Tensor< 1, dim > > &gradients) const; + + /** Return the gradient of the specified component of the function + at the given point. It is inefficient to use single points. If + you need more than one at a time, use the vector_value_list + function. For efficiency reasons, it is better if all the points + lie on the same cell. This is not mandatory, however it does + speed things up. */ + virtual Tensor<1,dim> gradient(const Point< dim > &p, + const unsigned int component = 0)const; + + /** Return the gradient of all components of the function at all + the given points. This is rather efficient if all the points + lie on the same cell. If this is not the case, things may slow + down a bit. */ + virtual void + vector_gradient_list (const std::vector< Point< dim > > &p, + std::vector< + std::vector< Tensor< 1, dim > > > &gradients) const; + + /** Return the gradient of the specified component of the function + at all the given points. This is rather efficient if all the + points lie on the same cell. If this is not the case, things + may slow down a bit. */ + virtual void + gradient_list (const std::vector< Point< dim > > &p, + std::vector< Tensor< 1, dim > > &gradients, + const unsigned int component=0) const; + + /** Create quadrature rules. This function groups the points into + blocks that live in the same cell, and fills up three vectors: + @p cells, @p qpoints and @p maps. The first is a list of the + cells that contain the points, the second is a list of + quadrature points matching each cell of the first list, and the + third contains the index of the given quadrature points, i.e., + @p points[maps[3][4]] ends up as the 5th quadrature point in the + 4th cell. This is where optimization would help. This function + returns the number of cells that contain the given set of + points. + */ + unsigned int compute_point_locations(const std::vector > &points, + std::vector &cells, + std::vector > > &qpoints, + std::vector > &maps) const; + + private: + /** Pointer to the dof handler. */ + SmartPointer dh; + + /** A reference to the actual data vector. */ + const VECTOR & data_vector; + + /** A reference to the mapping being used. */ + const Mapping & mapping; + + /** The current cell in which we are evaluating*/ + mutable typename DH::active_cell_iterator cell; + + /** Store the number of components of this function. */ + const unsigned int n_components; + }; +} +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/numerics/fe_field_function.templates.h b/deal.II/deal.II/include/numerics/fe_field_function.templates.h new file mode 100644 index 0000000000..9f07fb7f55 --- /dev/null +++ b/deal.II/deal.II/include/numerics/fe_field_function.templates.h @@ -0,0 +1,327 @@ +//--------------------------------------------------------------------------- +// $Id: function_parser.h 14594 2007-03-22 20:17:41Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2007 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions { + + template + FEFieldFunction::FEFieldFunction (const DH &mydh, + const VECTOR &myv, + const Mapping &mymapping) : + Function(mydh.get_fe().n_components()), + dh(&mydh, "FEFieldFunction"), + data_vector(myv), + mapping(mymapping), + n_components(mydh.get_fe().n_components()) + { + cell = dh->begin_active(); + } + + template + void FEFieldFunction::set_active_cell(typename DH::active_cell_iterator &newcell) { + cell = newcell; + } + + template + void FEFieldFunction::vector_value (const Point &p, + Vector &values) const + { + Assert (values.size() == n_components, + ExcDimensionMismatch(values.size(), n_components)); + 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)) { + std::pair > my_pair + = GridTools::find_active_cell_around_point (mapping, *dh, p); + cell = my_pair.first; + qp = my_pair.second; + } + + // Now we can find out about the point + Quadrature quad(qp); + FEValues fe_v(mapping, dh->get_fe(), quad, + update_values); + fe_v.reinit(cell); + std::vector< Vector > vvalues (1, values); + fe_v.get_function_values(data_vector, vvalues); + values = vvalues[0]; + } + + template + double FEFieldFunction::value + (const Point &p, unsigned int comp) const + { + Vector values(n_components); + vector_value(p, values); + return values(comp); + } + + + template + void FEFieldFunction::vector_gradient + (const Point &p, + std::vector > &gradients) const + { + Assert (gradients.size() == n_components, + ExcDimensionMismatch(gradients.size(), n_components)); + 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)) { + std::pair > my_pair + = GridTools::find_active_cell_around_point (mapping, *dh, p); + cell = my_pair.first; + qp = my_pair.second; + } + + // Now we can find out about the point + Quadrature quad(qp); + FEValues fe_v(mapping, dh->get_fe(), quad, + update_gradients); + fe_v.reinit(cell); + std::vector< std::vector > > vgrads + (1, std::vector >(n_components) ); + fe_v.get_function_grads(data_vector, vgrads); + gradients = vgrads[0]; + } + + template + Tensor<1,dim> FEFieldFunction::gradient + (const Point &p, unsigned int comp) const + { + std::vector > grads(n_components); + vector_gradient(p, grads); + return grads[comp]; + } + + // Now the list versions + // ============================== + + template + void FEFieldFunction::vector_value_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_values); + fe_v.reinit(cells[i]); + std::vector< Vector > vvalues (nq, Vector(n_components)); + fe_v.get_function_values(data_vector, vvalues); + for(unsigned int q=0; q + void FEFieldFunction::value_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_value_list(points, vvalues); + for(unsigned int q=0; q + void FEFieldFunction::vector_gradient_list (const std::vector > & points, + std::vector< + std::vector< Tensor<1,dim> > > &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_gradients); + fe_v.reinit(cells[i]); + std::vector< std::vector > > vgrads (nq, std::vector >(n_components)); + fe_v.get_function_grads(data_vector, vgrads); + for(unsigned int q=0; q + void FEFieldFunction::gradient_list (const std::vector > &points, + std::vector< Tensor<1,dim> > &values, + const unsigned int component) const + { + Assert(points.size() == values.size(), + ExcDimensionMismatch(points.size(), values.size())); + std::vector< std::vector > > vvalues(points.size(), std::vector >(n_components)); + vector_gradient_list(points, vvalues); + for(unsigned int q=0; q + unsigned int FEFieldFunction:: + compute_point_locations(const std::vector > &points, + std::vector &cells, + std::vector > > &qpoints, + std::vector > &maps) const { + // How many points are here? + unsigned int np = points.size(); + + // Reset output maps. + cells.clear(); + qpoints.clear(); + maps.clear(); + + // Now the easy case. + if(np==0) return 0; + + // Keep track of the points we found + std::vector point_flags(np, false); + + // Set this to true untill all points have been classified + bool left_over = true; + + // Current quadrature point + Point qp = mapping.transform_real_to_unit_cell(cell, points[0]); + + // Check if we already have a valid cell for the first point + if(!GeometryInfo::is_inside_unit_cell(qp)) { + std::pair > my_pair = GridTools::find_active_cell_around_point + (mapping, *dh, points[0]); + cell = my_pair.first; + qp = my_pair.second; + point_flags[0] = true; + } + + // Put in the first point. + cells.push_back(cell); + qpoints.push_back(std::vector >(1, qp)); + maps.push_back(std::vector (1, 0)); + + // Check if we need to do anything else + if(points.size() > 1) { + left_over = true; + } else { + left_over = false; + } + + // This is the first index of a non processed point + unsigned int first_outside = 1; + + // And this is the index of the current cell + unsigned int c = 0; + + while(left_over == true) { + // Assume this is the last one + left_over = false; + Assert(first_outside < np, + ExcIndexRange(first_outside, 0, np)); + + // If we found one in this cell, keep looking in the same cell + for(unsigned int p=first_outside; p qpoint = mapping.transform_real_to_unit_cell(cell, points[p]); + if(GeometryInfo::is_inside_unit_cell(qpoint)) { + point_flags[p] = true; + qpoints[c].push_back(qpoint); + maps[c].push_back(p); + } else { + // Set things up for next round + if(left_over == false) first_outside = p; + left_over = true; + } + } + // If we got here and there is no left over, we are done. Else we + // need to find the next cell + if(left_over == true) { + std::pair > my_pair + = GridTools::find_active_cell_around_point (mapping, *dh, points[first_outside]); + cells.push_back(my_pair.first); + qpoints.push_back(std::vector >(1, my_pair.second)); + maps.push_back(std::vector(1, first_outside)); + c++; + point_flags[first_outside] = true; + // And check if we can exit the loop now + if (first_outside == np-1) left_over = false; + } + } + + // Augment of one the number of cells + ++c; + // Debug Checking + Assert(c == cells.size(), ExcInternalError()); + + Assert(c == maps.size(), + ExcDimensionMismatch(c, maps.size())); + + Assert(c == qpoints.size(), + ExcDimensionMismatch(c, qpoints.size())); + +#ifdef DEBUG + unsigned int qps = 0; + // The number of points in all the cells must be the same as the + // number of points we started off from. + for(unsigned int n=0; n +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions { + + template class FEFieldFunction, + Vector >; + + template class FEFieldFunction, + BlockVector >; + + template class FEFieldFunction, + Vector >; + + template class FEFieldFunction, + BlockVector >; + +#ifdef DEAL_II_USE_PETSC + + template class FEFieldFunction, + PETScWrappers::Vector >; + + template class FEFieldFunction, + PETScWrappers::BlockVector >; + + template class FEFieldFunction, + PETScWrappers::Vector >; + + template class FEFieldFunction, + PETScWrappers::BlockVector >; + +#endif + +// template class FEFieldFunction, +// Vector >; + +// template class FEFieldFunction, +// BlockVector >; + +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/doc/authors.html b/deal.II/doc/authors.html index e4a40a724a..cf24a1e525 100644 --- a/deal.II/doc/authors.html +++ b/deal.II/doc/authors.html @@ -69,7 +69,7 @@ alphabetical order): generalization of FilteredMatrix;   integration of the function parser library;   cubit journal file to export to ucd mesh format;   - FEFunction and ParsedFunction classes;   + FEFieldFunction and ParsedFunction classes;   random bug fixes and enhancements.
  • Oliver Kayser-Herold: @@ -157,4 +157,4 @@ alphabetical order): - \ No newline at end of file + diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index 7c15ac8058..8de7544008 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -481,6 +481,13 @@ inconvenience this causes.
      +

    1. New: There is a new class: + Functions::FEFieldFunction which is a Function + interface to a finite element solution. +
      + (Luca Heltai 2007/08/29) + +

    2. Improved: FunctionDerivative is now derived from AutoDerivativeFunction and implements gradients as well, giving you automatic second derivatives of a function. @@ -1115,6 +1122,12 @@ inconvenience this causes.

      deal.II

        +

      1. New: There is a new class: + Functions::ParsedFunction which is friendly + wrapper to the FunctionParser class. +
        + (Luca Heltai 2007/08/29) +

      2. Fixed: the function DataOut::build_patches had a quadratic algorithm when generatic cell-data (as opposed -- 2.39.5