From 2faef03fddea37e96ce4feafe4a3b5bc6995d59e Mon Sep 17 00:00:00 2001 From: kanschat Date: Fri, 20 Oct 2006 15:54:12 +0000 Subject: [PATCH] reorganize FETools::get_fe_by_name and allow for adding custom elements git-svn-id: https://svn.dealii.org/trunk@14034 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_base.h | 7 + deal.II/deal.II/include/fe/fe_tools.h | 212 +++++++---- deal.II/deal.II/source/fe/fe_tools.cc | 485 ++++++++++++-------------- 3 files changed, 358 insertions(+), 346 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index c35f7edcaf..4a40d6a788 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -280,6 +280,13 @@ class FiniteElementData */ H2 = 0x0e }; + + /** + * The dimension of the finite + * element, which is the template + * parameter dim + */ + static const unsigned int dimension = dim; /** * Number of degrees of freedom on diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index 884a47c619..e8d54b79fe 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -21,6 +21,8 @@ #include +#include + template class FullMatrix; template class Vector; template class Quadrature; @@ -53,11 +55,56 @@ class ConstraintMatrix; * $id-I_h$ that is needed for evaluating $(id-I_h)z$ for e.g. the * dual solution $z$. * - * @author Ralf Hartmann, 2000; Wolfgang Bangerth, 2003, 2005, Guido Kanschat, 2000, 2004 + * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat; + * 2000, 2003, 2004, 2005, 2006 */ class FETools { public: + /** + * The base for factory objects + * creating finite elements of a + * given degree. + * + * @author Guido Kanschat, 2006 + */ + template + class FEFactoryBase + { + public: + /** + * Create a FiniteElement and + * return a pointer to it. + */ + virtual FiniteElement* + get (const unsigned int degree) const = 0; + /** + * Virtual destructor doing + * nothing but making the + * compiler happy. + */ + virtual ~FEFactoryBase(); + }; + + /** + * The base for factory objects + * creating finite elements of a + * given degree. + * + * @author Guido Kanschat, 2006 + */ + template + class FEFactory : public FEFactoryBase + { + public: + /** + * Create a FiniteElement and + * return a pointer to it. + */ + virtual FiniteElement* + get (const unsigned int degree) const; + }; + /** * @warning In most cases, you * will probably want to use @@ -909,54 +956,39 @@ class FETools std::vector &l2h); /** - * Given a name in the form which + * Parse the name of a finite + * element and generate a finite + * element object accordingly. + * + * The name must be in the form which * is returned by the * @p FiniteElement::get_name - * function, regenerate such a - * finite element. + * function, where a few + * modifications are allowed: + * + *
  • Dimension template + * parameters <2> etc. can + * be omitted. Alternatively, the + * explicit number can be + * replaced by dim or + * d. If a number is + * given, it must match + * the template parameter of this + * function. * - * This function is useful to - * convert the name given in an - * input file to an actual finite - * element, without having to - * parse the name yourself. + *
  • The powers used for + * FESystem may either be numbers + * or can be + * replaced by dim or + * d. + *
* - * Note that the given name must - * match exactly what one would - * get from the finite element to - * be created, since otherwise - * the parsing would fail. If no - * finite element can be + * If no finite element can be * reconstructed from this * string, an exception of type * @p FETools::ExcInvalidFEName * is thrown. * - * There is one exception, - * however, where the names must - * not match exactly: while the - * finite elements write the - * space dimension in the form of - * a template argument after the - * name of the class (for example - * FE_Q<2>, you can omit the - * dimension argument altogether, - * or replace it with the string - * @. The reason - * is that the dimension argument - * may be cumbersome if the name - * of a finite element is given - * in an input file that may be - * used to control operation of - * the program in different space - * dimensions. Running the - * program in another space - * dimension would then require - * changing the input file as - * well. With above exception, - * there is a canonical spelling - * that doesn't require this. - * * The function returns a pointer * to a newly create finite * element. It is in the callers @@ -976,21 +1008,56 @@ class FETools FiniteElement * get_fe_from_name (const std::string &name); - /** - * Exception + * Adds the name of a finite + * element to be used by + * get_fe_from_name(). + * + * It is safe to use either the + * class name explicitly or to + * use the result of + * FiniteElement::get_name, since + * everything after the first + * non-name character will be + * chopped off. + */ + template + static void add_fe_name(const std::string& name, + boost::shared_ptr > factory); + + /** + * The string used for + * get_fe_from_name() cannot be + * translated to a finite + * element. + * + * Either the string is badly + * formatted or you are using a + * custom element that must be + * added using add_fe_name() + * first. */ DeclException1 (ExcInvalidFEName, std::string, - << "Can't re-generate a finite element from the string <" - << arg1 << ">."); + << "Can't re-generate a finite element from the string '" + << arg1 << "'."); + + /** + * Parsing a finite element name, + * an unexpected character showed up. + */ + DeclException1 (ExcInvalidFECharacter, + std::string, + << "Unexpected character at beginning of '"<< arg1 << "'"); + /** * Exception */ DeclException0 (ExcInvalidFE); /** - * Exception + * The finite element must be + * @ref GlossPrimitive "primitive". */ DeclException0 (ExcFEMustBePrimitive); /** @@ -1008,7 +1075,7 @@ class FETools << "constraints. Use the respective function with " << "additional ConstraintMatrix argument(s), instead."); /** - * Exception + * You need at least two grid levels. */ DeclException0 (ExcGridNotRefinedAtLeastOnce); /** @@ -1019,36 +1086,6 @@ class FETools << "This is a " << arg1 << "x" << arg2 << " matrix, " << "but should be a " << arg3 << "x" << arg4 << " matrix."); - private: - /** - * Return a finite element that - * is created using the - * characters of the input - * parameters. The second part of - * the return value indicates how - * many characters have been used - * up in the creation of the - * finite element, so that the - * calling site can continue - * parsing finite element lists - * (for example for - * FESystem objects) at the - * position after which the - * present element's name ends. - * - * If no finite element could be - * created from the string at the - * beginning of the given string, - * then an exception is thrown, - * just as for the - * get_fe_from_name() - * function. - */ - template - static - std::pair *, unsigned int> - get_fe_from_name_aux (const std::string &name); - /** * Exception thrown if an * embedding matrix was computed @@ -1056,8 +1093,31 @@ class FETools */ DeclException1(ExcLeastSquaresError, double, << "Least squares fit leaves a gap of " << arg1); + private: + /** + * Return a finite element that + * is created using the beginning + * of name and eat away + * the part of name + * defining this element. + */ + template + static + FiniteElement * + get_fe_from_name_aux (std::string &name); + }; + +template +FiniteElement* +FETools::FEFactory::get (const unsigned int degree) const +{ + return new FE(degree); +} + + + /*@}*/ /*---------------------------- fe_tools.h ---------------------------*/ diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 2baa4c5ae6..d423671632 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -11,6 +11,7 @@ // //--------------------------------------------------------------------------- +#include #include #include @@ -29,6 +30,7 @@ #include #include #include +#include #include #include #include @@ -42,6 +44,16 @@ #include +// This is the map used by FETools::get_fe_from_name and +// FETools::add_fe_name. Since FEFactoryBase has a template parameter +// dim, it could not be a member variable of FETools. On the other +// hand, it is only accessed by functions in this file, so it is safe +// to make it a static variable here. It must be static so that we can +// link several dimensions together. +static std::map > > +fe_name_map; + namespace { @@ -155,6 +167,10 @@ namespace } +template +FETools::FEFactoryBase::~FEFactoryBase() +{} + template void FETools::compute_component_wise( @@ -1369,283 +1385,188 @@ void FETools::extrapolate(const DoFHandler &dof1, } - template -FiniteElement * -FETools::get_fe_from_name (const std::string &name) +void +FETools::add_fe_name(const std::string& parameter_name, + boost::shared_ptr > factory) { - // get the finite element that - // would be created from the given - // string at the first position, as - // well as how many characters of - // the string were eaten - const std::pair*, unsigned int> - tmp = get_fe_from_name_aux (name); - // make sure that we took all - // characters in the name, - // i.e. that there is not some junk - // left over at the end of which we - // didn't know what to do - // with. make sure we don't create - // a memory leak here - if (tmp.second != name.size()) - { - delete tmp.first; - AssertThrow (false, ExcInvalidFEName (name)); - } - - // otherwise, return the just - // created pointer - return tmp.first; + // Erase everything after the + // actual class name + std::string name = parameter_name; + unsigned int name_end = + name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); + if (name_end < name.size()) + name.erase(name_end); + // Insert the normalized name into + // the map + fe_name_map.insert(std::make_pair(name, factory)); } - template -std::pair *, unsigned int> -FETools::get_fe_from_name_aux (const std::string &name) -{ - // so, let's see what's at position - // 0 of this string, and create a - // respective finite element - // - // start with the longest names, to - // make sure we don't match FE_Q - // when it's actually a - // FE_Q_Hierarchic - if (Utilities::match_at_string_start (name, std::string("FE_Q_Hierarchical"))) - { - unsigned int position = std::string("FE_Q_Hierarchical").size(); - // as described in the - // documentation, at this point - // we may have either a) no - // dimension specification, b) - // the correct dimension (like - // <2>), or c) a generic - // . check how many of - // these characters match, and - // advance the cursor by the - // respective amount - position += match_dimension (name, position); - - // make sure the next character - // is an opening parenthesis - AssertThrow (name[position] == '(', ExcInvalidFEName(name)); - ++position; - // next thing is to parse the - // degree of the finite element - const std::pair tmp - = Utilities::get_integer_at_position (name, position); - AssertThrow (tmp.first>=0, ExcInvalidFEName(name)); - position += tmp.second; - - // make sure the next character - // is an closing parenthesis - AssertThrow (name[position] == ')', ExcInvalidFEName(name)); - ++position; - - // ok, everything seems - // good. so create finite - // element and return position - // count - return std::make_pair (static_cast*> - (new FE_Q_Hierarchical(tmp.first)), - position); - } - // check other possibilities in - // exactly the same way - else if (Utilities::match_at_string_start (name, std::string("FE_RaviartThomas"))) - { - unsigned int position = std::string("FE_RaviartThomas").size(); - position += match_dimension (name, position); - AssertThrow (name[position] == '(', ExcInvalidFEName(name)); - ++position; - - const std::pair tmp - = Utilities::get_integer_at_position (name, position); - - AssertThrow (tmp.first>=0, ExcInvalidFEName(name)); - position += tmp.second; - AssertThrow (name[position] == ')', ExcInvalidFEName(name)); - ++position; - return std::make_pair (static_cast*> - (new FE_RaviartThomas(tmp.first)), - position); - } - // check other possibilities in - // exactly the same way - else if (Utilities::match_at_string_start (name, std::string("FE_RaviartThomasNodal"))) +FiniteElement * +FETools::get_fe_from_name (const std::string ¶meter_name) +{ + // First, make sure that the + // standard map from names to + // element factories has been + // filled by checking for one of + // the elements. +//TODO: Make this threadsafe + if (fe_name_map.find(std::string("FE_Q_Hierarchical")) + == fe_name_map.end()) { - unsigned int position = std::string("FE_RaviartThomasNodal").size(); - position += match_dimension (name, position); - AssertThrow (name[position] == '(', ExcInvalidFEName(name)); - ++position; - - const std::pair tmp - = Utilities::get_integer_at_position (name, position); - - AssertThrow (tmp.first>=0, ExcInvalidFEName(name)); - position += tmp.second; - AssertThrow (name[position] == ')', ExcInvalidFEName(name)); - ++position; - return std::make_pair (static_cast*> - (new FE_RaviartThomasNodal(tmp.first)), - position); + add_fe_name (std::string("FE_Q_Hierarchical"), + boost::shared_ptr > + (new FEFactory >())); + add_fe_name (std::string("FE_ABF"), + boost::shared_ptr > + (new FEFactory >())); + add_fe_name (std::string("FE_RaviartThomas"), + boost::shared_ptr > + (new FEFactory >())); + add_fe_name (std::string("FE_RaviartThomasNodal"), + boost::shared_ptr > + (new FEFactory >())); + add_fe_name (std::string("FE_Nedelec"), + boost::shared_ptr > + (new FEFactory >())); + add_fe_name (std::string("FE_DGPNonparametric"), + boost::shared_ptr > + (new FEFactory >())); + add_fe_name (std::string("FE_DGP"), + boost::shared_ptr > + (new FEFactory >())); + add_fe_name (std::string("FE_DGPMonomial"), + boost::shared_ptr > + (new FEFactory >())); + add_fe_name (std::string("FE_DGQ"), + boost::shared_ptr > + (new FEFactory >())); + add_fe_name (std::string("FE_Q<2>(4)"), + boost::shared_ptr > + (new FEFactory >())); } - else if (Utilities::match_at_string_start (name, std::string("FE_Nedelec"))) + + // Create a version of the name + // string where all template + // parameters are eliminated. + std::string name = parameter_name; + for (unsigned int pos1 = name.find('<'); + pos1 < name.size(); + pos1 = name.find('<')) { - unsigned int position = std::string("FE_Nedelec").size(); - position += match_dimension (name, position); - AssertThrow (name[position] == '(', ExcInvalidFEName(name)); - ++position; - - const std::pair tmp - = Utilities::get_integer_at_position (name, position); - AssertThrow (tmp.first>=0, ExcInvalidFEName(name)); - position += tmp.second; - AssertThrow (name[position] == ')', ExcInvalidFEName(name)); - ++position; - return std::make_pair (static_cast*> - (new FE_Nedelec(tmp.first)), - position); - } - else if (Utilities::match_at_string_start (name, std::string("FE_DGPNonparametric"))) - { - unsigned int position = std::string("FE_DGPNonparametric").size(); - position += match_dimension (name, position); - AssertThrow (name[position] == '(', ExcInvalidFEName(name)); - ++position; - - const std::pair tmp - = Utilities::get_integer_at_position (name, position); + const unsigned int pos2 = name.find('>'); + // If there is only a single + // character between those two, + // it should be 'd' or the number + // representing the dimension. + if (pos2-pos1 == 2) + { + const char dimchar = '0' + dim; +//TODO: Write a more explicit exception here! + if (name.at(pos1+1) != 'd') + Assert (name.at(pos1+1) == dimchar, + ExcInvalidFEName(name)); + } else + Assert(pos2-pos1 == 4, ExcInvalidFEName(name)); - AssertThrow (tmp.first>=0, ExcInvalidFEName(name)); - position += tmp.second; - AssertThrow (name[position] == ')', ExcInvalidFEName(name)); - ++position; - return std::make_pair (static_cast*> - (new FE_DGPNonparametric(tmp.first)), - position); + // If pos1==pos2, then we are + // probably at the end of the + // string + if (pos2 != pos1) + name.erase(pos1, pos2-pos1+1); } - else if (Utilities::match_at_string_start (name, std::string("FE_DGP"))) + // Replace all occurences of "^dim" + // by "^d" to be handled by the + // next loop + for (unsigned int pos = name.find("^dim"); + pos < name.size(); + pos = name.find("^dim")) { - unsigned int position = std::string("FE_DGP").size(); - position += match_dimension (name, position); - AssertThrow (name[position] == '(', ExcInvalidFEName(name)); - ++position; - - const std::pair tmp - = Utilities::get_integer_at_position (name, position); - - AssertThrow (tmp.first>=0, ExcInvalidFEName(name)); - position += tmp.second; - AssertThrow (name[position] == ')', ExcInvalidFEName(name)); - ++position; - return std::make_pair (static_cast*>(new FE_DGP(tmp.first)), - position); + name.erase(pos+2, 2); } - else if (Utilities::match_at_string_start (name, std::string("FE_DGQ"))) + + // Replace all occurences of "^d" + // by using the actual dimension + for (unsigned int pos = name.find("^d"); + pos < name.size(); + pos = name.find("^d")) { - unsigned int position = std::string("FE_DGQ").size(); - position += match_dimension (name, position); - AssertThrow (name[position] == '(', ExcInvalidFEName(name)); - ++position; - - const std::pair tmp - = Utilities::get_integer_at_position (name, position); - - AssertThrow (tmp.first>=0, ExcInvalidFEName(name)); - position += tmp.second; - AssertThrow (name[position] == ')', ExcInvalidFEName(name)); - ++position; - return std::make_pair (static_cast*> (new FE_DGQ(tmp.first)), - position); + name.at(pos+1) = '0' + dim; } - else if (Utilities::match_at_string_start (name, std::string("FE_Q"))) - { - unsigned int position = std::string("FE_Q").size(); - position += match_dimension (name, position); - AssertThrow (name[position] == '(', ExcInvalidFEName(name)); - ++position; + FiniteElement* fe = get_fe_from_name_aux (name); + // Make sure the auxiliary function + // ate up all characters of the name. + AssertThrow (name.size() == 0, ExcInvalidFEName(parameter_name)); + return fe; +} - const std::pair tmp - = Utilities::get_integer_at_position (name, position); - - AssertThrow (tmp.first>=0, ExcInvalidFEName(name)); - position += tmp.second; - AssertThrow (name[position] == ')', ExcInvalidFEName(name)); - ++position; - return std::make_pair (static_cast*>(new FE_Q(tmp.first)), - position); - } + + +template +FiniteElement* +FETools::get_fe_from_name_aux (std::string &name) +{ + // Extract the name of the + // finite element class, which only + // contains characters, numbers and + // underscores. + unsigned int name_end = + name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); + const std::string name_part(name, 0, name_end); + name.erase(0, name_part.size()); - else - // now things get a little more - // complicated: FESystem. it's - // more complicated, since we - // have to figure out what the - // base elements are. this can - // only be done recursively - if (Utilities::match_at_string_start (name, std::string("FESystem"))) - { - unsigned int position = std::string("FESystem").size(); - position += match_dimension (name, position); - - // FESystem puts the names of - // the basis elements into - // square brackets - AssertThrow (name[position] == '[', ExcInvalidFEName(name)); - ++position; - - // next we have to get at the - // base elements. start with - // the first. wrap the whole - // block into try-catch to - // make sure we destroy the - // pointers we got from - // recursive calls if one of - // these calls should throw - // an exception - std::vector*> base_fes; - std::vector base_multiplicities; - try - { - do - { - // first get the - // element at this - // position of the - // string, i.e. of - // the substring - // ranging from the - // present position - // to the end - const std::pair *, unsigned int> tmp_x - = get_fe_from_name_aux (std::string(name.begin()+position, - name.end())); - base_fes.push_back (tmp_x.first); - position += tmp_x.second; - - // next check whether - // FESystem placed a - // multiplicity after - // the element name - if (name[position] == '^') - { - // yes. this is - // the case. move - // the cursor - // beyond the '^' - // and read this - // multiplicity - ++position; - - const std::pair tmp - = Utilities::get_integer_at_position (name, position); - - AssertThrow (tmp.first>=0, ExcInvalidFEName(name)); - position += tmp.second; + // now things get a little more + // complicated: FESystem. it's + // more complicated, since we + // have to figure out what the + // base elements are. this can + // only be done recursively + if (name_part == "FESystem") + { + // next we have to get at the + // base elements. start with + // the first. wrap the whole + // block into try-catch to + // make sure we destroy the + // pointers we got from + // recursive calls if one of + // these calls should throw + // an exception + std::vector*> base_fes; + std::vector base_multiplicities; + try + { + // Now, just the [...] + // part should be left. + AssertThrow(name.at(0) == '[', ExcInvalidFECharacter(name)); + do + { + // Erase the + // leading '[' or '-' + name.erase(0,1); + // Now, the name of the + // first base element is + // first... Let's get it + base_fes.push_back (get_fe_from_name_aux (name)); + // next check whether + // FESystem placed a + // multiplicity after + // the element name + if (name[0] == '^') + { + // yes. Delete the '^' + // and read this + // multiplicity + name.erase(0,1); + + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + name.erase(0, tmp.second); + // add to length, + // including the '^' base_multiplicities.push_back (tmp.first); } else @@ -1664,16 +1585,16 @@ FETools::get_fe_from_name_aux (const std::string &name) // so loop while the // next character is // '-' - } - while (name[position++] == '-'); + } + while (name[0] == '-'); // so we got to the end // of the '-' separated // list. make sure that // we actually had a ']' // there - AssertThrow (name[position-1] == ']', ExcInvalidFEName(name)); - + AssertThrow (name.at(0) == ']', ExcInvalidFECharacter(name)); + name.erase(0,1); // just one more sanity check Assert ((base_fes.size() == base_multiplicities.size()) && @@ -1715,7 +1636,7 @@ FETools::get_fe_from_name_aux (const std::string &name) } default: - Assert (false, ExcNotImplemented()); + AssertThrow (false, ExcNotImplemented()); } // now we don't need the @@ -1724,8 +1645,12 @@ FETools::get_fe_from_name_aux (const std::string &name) for (unsigned int i=0; i > >::const_iterator + entry = fe_name_map.find(name_part); + AssertThrow (entry != fe_name_map.end(), ExcInvalidFEName(name)); + // Now, just the (degree) + // part should be left. + AssertThrow(name.at(0) == '(', ExcInvalidFECharacter(name)); + name.erase(0,1); + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + name.erase(0, tmp.second+1); + return entry->second->get(tmp.first); + } // hm, if we have come thus far, we @@ -1760,7 +1700,7 @@ FETools::get_fe_from_name_aux (const std::string &name) // make some compilers happy that // do not realize that we can't get // here after throwing - return std::pair *, unsigned int> (0,0); + return 0; } @@ -1829,6 +1769,8 @@ compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, /*-------------- Explicit Instantiations -------------------------------*/ +template class FETools::FEFactoryBase; + template void FETools::compute_node_matrix( FullMatrix&, @@ -2085,10 +2027,13 @@ void FETools::extrapolate const DoFHandler &, const ConstraintMatrix &, Vector &); -template -FiniteElement * +template FiniteElement * FETools::get_fe_from_name (const std::string &); - +template void +FETools::add_fe_name( + const std::string& name, + boost::shared_ptr > factory); + template void FETools:: -- 2.39.5