The deal.II library has been placed under an Open
Source license, in the sense advocated by the Open Source
- Initiative. You are thus free to copy and use it, and you have
- free access to all source codes. However, deal.II
- is not in the public domain, it is property of and copyrighted by
- the deal.II authors, and there are restrictions
- on its use. We will give some hints on license issues first. The
- legally binding license is the Q Public
- License (QPL) at the bottom of this page. Furthermore, deal.II
- distributions come bundled with some third party software for which
- the deal.II authors do not own the copyright. The licenses of these
- parts are discussed below.
+ href="http://www.opensource.org" target="_top">Open Source
+ Initiative. You are thus free to copy and use it, and you have
+ free access to all source code.
+
+
+ However, deal.II is not in the public domain, it
+ is property of and copyrighted by the
+ deal.II authors, and
+ there are restrictions on its use:
+ The legally binding license is the GNU Lesser
+ General Public License (LGPL) as published by the Free Software
+ Foundation; either version 2.1 of the License, or (at your option)
+ any later version.
+
+
+
+ Furthermore, the deal.II distribution comes bundled with some third
+ party software for which the deal.II authors do not own the
+ copyright. The licenses of these parts are discussed
+ below.
+
In this document, as well as in the program files, the terms
@@ -40,247 +50,609 @@
Kanschat and a long list of contributors.
-
-
-
Restrictions on and commercial use of the library
-
-
-
- deal.II is distributed under the Open Source
- license called the Q Public License (or "QPL") listed
- in full below. Essentially, an interpretation of this license is
- that you can use the library free of charge for developing your own
- programs. However, if you plan to distribute these programs to others,
- certain restrictions are imposed by the QPL.
+
Contributions
-
- In non-legal terms, the license stipulates that if you want
- to give your program to someone else in whatever form (including as
- a binary executable), you have to make the entire source code of
- your own program and all modifications to the deal.II library
- available free of charge as well.
-
+ As a contributor to this project, you agree that any of your
+ contributions be licensed under the same terms and conditions as
+ the license of the deal.II project granted to you.
+
-
- We recognize that giving away the source code may be an
- unacceptable restriction in a commercial setting. As the owners of
- the sources, we are able to grant you a differently licensed version
- of the deal.II library that does not contain this
- requirement, and we are very much willing to negotiate the terms of
- such a license. Please contact the authors of the library if you are
- interested in commercial use of deal.II.
-
+ We, the deal.II authors, do not require
+ copyright assignments for contributions. This means that the
+ copyright for code contributions in the deal.II project is held by
+ its respective contributors who have each agreed to release their
+ contributed code under the terms of the LGPL v2.1 or later.
+
-
Referencing the library
+
Referencing the library
- In addition to the terms imposed by the QPL, we ask for the
- following courtesy:
-
+ In addition to the terms imposed by the LGPL v2.1 or later, we ask
+ for the following courtesy:
-
Every publication presenting numerical results obtained with
- the help of deal.II should state the name of
- the library and cite
- the canonical deal.II reference listed at the top of the publications page.
+
Every publication presenting numerical results obtained with the
+ help of deal.II should state the name of the
+ library and cite the canonical deal.II reference listed at the top of
+ the publications page.
- This is the usual, fair way of giving credit to contributors to a
- scientific result. In addition, it helps us justify our effort in
- developing deal.II as an academic undertaking.
-
+
+
+ This is the usual, fair way of giving credit to contributors to a
+ scientific result. In addition, it helps us justify our effort in
+ developing deal.II as an academic undertaking.
+
-
-
- Some of the files in the contrib/ subdirectory are by
- other people, and are distributed under separate licenses. These
- files are not part of the deal.II library, but
- they are part of our distribution. They are generally available for
- use in non-commercial academic settings. For commercial use and if
- you want to distribute your codes, please consult the licenses of
- these parts of the deal.II distribution. We made sure that those parts of the
- distribution for which the license is not entirely free for all uses can be
- disabled upon configuration of deal.II without
- unduly restricting the functionality of the rest of the
- library.
-
-
-
- In particular, these are libraries that we distribute in the
- contrib/ subdirectory, courtesy of their authors:
-
- The QPL was originally developed by
- Troll Tech for their
- Qt program library, but has since gained some support by other
- library authors as well.
- For more information on the QPL, please follow
- this link.
+ deal.II comes bundled with copies of
+ UMFPACK,
+ Threading Building Blocks,
+ BOOST and
+ FunctionParser.
+ These projects reside in the subdirectory "bundled/". Please note
+ that these projects are copyrighted by others than the deal.II
+ authors and are covered by different licenses. For details,
+ consult the respective webpage.
- THE Q PUBLIC LICENSE
- version 1.0
-
-
-
- Copyright (C) 1999 Troll Tech AS, Norway.
-
-
- Everyone is permitted to copy and
-
- distribute this license document.
-
-
-
- The intent of this license is to establish freedom to share and change the
- software regulated by this license under the open source model.
-
- This license applies to any software containing a notice placed by the
- copyright holder saying that it may be distributed under the terms of
- the Q Public License version 1.0. Such software is herein referred to as
- the Software. This license covers modification and distribution of the
- Software, use of third-party application programs based on the Software,
- and development of free software which uses the Software.
-
- Granted Rights
-
-
- 1. You are granted the non-exclusive rights set forth in this license
- provided you agree to and comply with any and all conditions in this
- license. Whole or partial distribution of the Software, or software
- items that link with the Software, in any form signifies acceptance of
- this license.
-
- 2. You may copy and distribute the Software in unmodified form provided
- that the entire package, including - but not restricted to - copyright,
- trademark notices and disclaimers, as released by the initial developer
- of the Software, is distributed.
-
- 3. You may make modifications to the Software and distribute your
- modifications, in a form that is separate from the Software, such as
- patches. The following restrictions apply to modifications:
-
-
- a. Modifications must not alter or remove any copyright notices in
- the Software.
-
- b. When modifications to the Software are released under this
- license, a non-exclusive royalty-free right is granted to the
- initial developer of the Software to distribute your modification
- in future versions of the Software provided such versions remain
- available under these terms in addition to any other license(s) of
- the initial developer.
-
-
- 4. You may distribute machine-executable forms of the Software or
- machine-executable forms of modified versions of the Software, provided
- that you meet these restrictions:
-
-
- a. You must include this license document in the distribution.
-
- b. You must ensure that all recipients of the machine-executable forms
- are also able to receive the complete machine-readable source code
- to the distributed Software, including all modifications, without
- any charge beyond the costs of data transfer, and place prominent
- notices in the distribution explaining this.
-
- c. You must ensure that all modifications included in the
- machine-executable forms are available under the terms of this
- license.
-
-
- 5. You may use the original or modified versions of the Software to
- compile, link and run application programs legally developed by you
- or by others.
-
- 6. You may develop application programs, reusable components and other
- software items that link with the original or modified versions of the
- Software. These items, when distributed, are subject to the following
- requirements:
-
-
- a. You must ensure that all recipients of machine-executable forms of
- these items are also able to receive and use the complete
- machine-readable source code to the items without any charge
- beyond the costs of data transfer.
-
- b. You must explicitly license all recipients of your items to use
- and re-distribute original and modified versions of the items in
- both machine-executable and source code forms. The recipients must
- be able to do so without any charges whatsoever, and they must be
- able to re-distribute to anyone they choose.
-
- c. If the items are not available to the general public, and the
- initial developer of the Software requests a copy of the items,
- then you must supply one.
-
-
- Limitations of Liability
-
- In no event shall the initial developers or copyright holders be liable
- for any damages whatsoever, including - but not restricted to - lost
- revenue or profits or other direct, indirect, special, incidental or
- consequential damages, even if they have been advised of the possibility
- of such damages, except to the extent invariable law, if any, provides
- otherwise.
-
- No Warranty
-
+ The libraries listed above are all open source and as such place few
+ restrictions on their use. However, if you want to distribute
+ software developed with deal.II (in source or binary form) and you
+ are using the packages above (with source code in
+ bundled/), then they may impose different terms. Please
+ consult the licenses of these packages for more information.
+
- The Software and this license document are provided AS IS with NO WARRANTY
- OF ANY KIND, INCLUDING THE WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS
- FOR A PARTICULAR PURPOSE.
-
- Choice of Law
-
+ Alternatively, the configuration process of deal.II allows you to
+ disable the use of all of these bundled libraries. In that case, you
+ can not use their functionality but you also are not restricted by
+ their license.
- This license is governed by the Laws of Germany. Disputes shall be settled
- by Amtsgericht Heidelberg.
+
+
Contact
+
+ For further questions regarding licensing and commercial use please
+ contact the deal.II maintainers directly via one of the following mail
+ addresses:
+
+
+ GNU LESSER GENERAL PUBLIC LICENSE
+ Version 2.1, February 1999
+
+ Copyright (C) 1991, 1999 Free Software Foundation, Inc.
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ [This is the first released version of the Lesser GPL. It also counts
+ as the successor of the GNU Library Public License, version 2, hence
+ the version number 2.1.]
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+ freedom to share and change it. By contrast, the GNU General Public
+ Licenses are intended to guarantee your freedom to share and change
+ free software--to make sure the software is free for all its users.
+
+ This license, the Lesser General Public License, applies to some
+ specially designated software packages--typically libraries--of the
+ Free Software Foundation and other authors who decide to use it. You
+ can use it too, but we suggest you first think carefully about whether
+ this license or the ordinary General Public License is the better
+ strategy to use in any particular case, based on the explanations below.
+
+ When we speak of free software, we are referring to freedom of use,
+ not price. Our General Public Licenses are designed to make sure that
+ you have the freedom to distribute copies of free software (and charge
+ for this service if you wish); that you receive source code or can get
+ it if you want it; that you can change the software and use pieces of
+ it in new free programs; and that you are informed that you can do
+ these things.
+
+ To protect your rights, we need to make restrictions that forbid
+ distributors to deny you these rights or to ask you to surrender these
+ rights. These restrictions translate to certain responsibilities for
+ you if you distribute copies of the library or if you modify it.
+
+ For example, if you distribute copies of the library, whether gratis
+ or for a fee, you must give the recipients all the rights that we gave
+ you. You must make sure that they, too, receive or can get the source
+ code. If you link other code with the library, you must provide
+ complete object files to the recipients, so that they can relink them
+ with the library after making changes to the library and recompiling
+ it. And you must show them these terms so they know their rights.
+
+ We protect your rights with a two-step method: (1) we copyright the
+ library, and (2) we offer you this license, which gives you legal
+ permission to copy, distribute and/or modify the library.
+
+ To protect each distributor, we want to make it very clear that
+ there is no warranty for the free library. Also, if the library is
+ modified by someone else and passed on, the recipients should know
+ that what they have is not the original version, so that the original
+ author's reputation will not be affected by problems that might be
+ introduced by others.
+
+ Finally, software patents pose a constant threat to the existence of
+ any free program. We wish to make sure that a company cannot
+ effectively restrict the users of a free program by obtaining a
+ restrictive license from a patent holder. Therefore, we insist that
+ any patent license obtained for a version of the library must be
+ consistent with the full freedom of use specified in this license.
+
+ Most GNU software, including some libraries, is covered by the
+ ordinary GNU General Public License. This license, the GNU Lesser
+ General Public License, applies to certain designated libraries, and
+ is quite different from the ordinary General Public License. We use
+ this license for certain libraries in order to permit linking those
+ libraries into non-free programs.
+
+ When a program is linked with a library, whether statically or using
+ a shared library, the combination of the two is legally speaking a
+ combined work, a derivative of the original library. The ordinary
+ General Public License therefore permits such linking only if the
+ entire combination fits its criteria of freedom. The Lesser General
+ Public License permits more lax criteria for linking other code with
+ the library.
+
+ We call this license the "Lesser" General Public License because it
+ does Less to protect the user's freedom than the ordinary General
+ Public License. It also provides other free software developers Less
+ of an advantage over competing non-free programs. These disadvantages
+ are the reason we use the ordinary General Public License for many
+ libraries. However, the Lesser license provides advantages in certain
+ special circumstances.
+
+ For example, on rare occasions, there may be a special need to
+ encourage the widest possible use of a certain library, so that it becomes
+ a de-facto standard. To achieve this, non-free programs must be
+ allowed to use the library. A more frequent case is that a free
+ library does the same job as widely used non-free libraries. In this
+ case, there is little to gain by limiting the free library to free
+ software only, so we use the Lesser General Public License.
+
+ In other cases, permission to use a particular library in non-free
+ programs enables a greater number of people to use a large body of
+ free software. For example, permission to use the GNU C Library in
+ non-free programs enables many more people to use the whole GNU
+ operating system, as well as its variant, the GNU/Linux operating
+ system.
+
+ Although the Lesser General Public License is Less protective of the
+ users' freedom, it does ensure that the user of a program that is
+ linked with the Library has the freedom and the wherewithal to run
+ that program using a modified version of the Library.
+
+ The precise terms and conditions for copying, distribution and
+ modification follow. Pay close attention to the difference between a
+ "work based on the library" and a "work that uses the library". The
+ former contains code derived from the library, whereas the latter must
+ be combined with the library in order to run.
+
+ GNU LESSER GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License Agreement applies to any software library or other
+ program which contains a notice placed by the copyright holder or
+ other authorized party saying it may be distributed under the terms of
+ this Lesser General Public License (also called "this License").
+ Each licensee is addressed as "you".
+
+ A "library" means a collection of software functions and/or data
+ prepared so as to be conveniently linked with application programs
+ (which use some of those functions and data) to form executables.
+
+ The "Library", below, refers to any such software library or work
+ which has been distributed under these terms. A "work based on the
+ Library" means either the Library or any derivative work under
+ copyright law: that is to say, a work containing the Library or a
+ portion of it, either verbatim or with modifications and/or translated
+ straightforwardly into another language. (Hereinafter, translation is
+ included without limitation in the term "modification".)
+
+ "Source code" for a work means the preferred form of the work for
+ making modifications to it. For a library, complete source code means
+ all the source code for all modules it contains, plus any associated
+ interface definition files, plus the scripts used to control compilation
+ and installation of the library.
+
+ Activities other than copying, distribution and modification are not
+ covered by this License; they are outside its scope. The act of
+ running a program using the Library is not restricted, and output from
+ such a program is covered only if its contents constitute a work based
+ on the Library (independent of the use of the Library in a tool for
+ writing it). Whether that is true depends on what the Library does
+ and what the program that uses the Library does.
+
+ 1. You may copy and distribute verbatim copies of the Library's
+ complete source code as you receive it, in any medium, provided that
+ you conspicuously and appropriately publish on each copy an
+ appropriate copyright notice and disclaimer of warranty; keep intact
+ all the notices that refer to this License and to the absence of any
+ warranty; and distribute a copy of this License along with the
+ Library.
+
+ You may charge a fee for the physical act of transferring a copy,
+ and you may at your option offer warranty protection in exchange for a
+ fee.
+
+ 2. You may modify your copy or copies of the Library or any portion
+ of it, thus forming a work based on the Library, and copy and
+ distribute such modifications or work under the terms of Section 1
+ above, provided that you also meet all of these conditions:
+
+ a) The modified work must itself be a software library.
+
+ b) You must cause the files modified to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ c) You must cause the whole of the work to be licensed at no
+ charge to all third parties under the terms of this License.
+
+ d) If a facility in the modified Library refers to a function or a
+ table of data to be supplied by an application program that uses
+ the facility, other than as an argument passed when the facility
+ is invoked, then you must make a good faith effort to ensure that,
+ in the event an application does not supply such function or
+ table, the facility still operates, and performs whatever part of
+ its purpose remains meaningful.
+
+ (For example, a function in a library to compute square roots has
+ a purpose that is entirely well-defined independent of the
+ application. Therefore, Subsection 2d requires that any
+ application-supplied function or table used by this function must
+ be optional: if the application does not supply it, the square
+ root function must still compute square roots.)
+
+ These requirements apply to the modified work as a whole. If
+ identifiable sections of that work are not derived from the Library,
+ and can be reasonably considered independent and separate works in
+ themselves, then this License, and its terms, do not apply to those
+ sections when you distribute them as separate works. But when you
+ distribute the same sections as part of a whole which is a work based
+ on the Library, the distribution of the whole must be on the terms of
+ this License, whose permissions for other licensees extend to the
+ entire whole, and thus to each and every part regardless of who wrote
+ it.
+
+ Thus, it is not the intent of this section to claim rights or contest
+ your rights to work written entirely by you; rather, the intent is to
+ exercise the right to control the distribution of derivative or
+ collective works based on the Library.
+
+ In addition, mere aggregation of another work not based on the Library
+ with the Library (or with a work based on the Library) on a volume of
+ a storage or distribution medium does not bring the other work under
+ the scope of this License.
+
+ 3. You may opt to apply the terms of the ordinary GNU General Public
+ License instead of this License to a given copy of the Library. To do
+ this, you must alter all the notices that refer to this License, so
+ that they refer to the ordinary GNU General Public License, version 2,
+ instead of to this License. (If a newer version than version 2 of the
+ ordinary GNU General Public License has appeared, then you can specify
+ that version instead if you wish.) Do not make any other change in
+ these notices.
+
+ Once this change is made in a given copy, it is irreversible for
+ that copy, so the ordinary GNU General Public License applies to all
+ subsequent copies and derivative works made from that copy.
+
+ This option is useful when you wish to copy part of the code of
+ the Library into a program that is not a library.
+
+ 4. You may copy and distribute the Library (or a portion or
+ derivative of it, under Section 2) in object code or executable form
+ under the terms of Sections 1 and 2 above provided that you accompany
+ it with the complete corresponding machine-readable source code, which
+ must be distributed under the terms of Sections 1 and 2 above on a
+ medium customarily used for software interchange.
+
+ If distribution of object code is made by offering access to copy
+ from a designated place, then offering equivalent access to copy the
+ source code from the same place satisfies the requirement to
+ distribute the source code, even though third parties are not
+ compelled to copy the source along with the object code.
+
+ 5. A program that contains no derivative of any portion of the
+ Library, but is designed to work with the Library by being compiled or
+ linked with it, is called a "work that uses the Library". Such a
+ work, in isolation, is not a derivative work of the Library, and
+ therefore falls outside the scope of this License.
+
+ However, linking a "work that uses the Library" with the Library
+ creates an executable that is a derivative of the Library (because it
+ contains portions of the Library), rather than a "work that uses the
+ library". The executable is therefore covered by this License.
+ Section 6 states terms for distribution of such executables.
+
+ When a "work that uses the Library" uses material from a header file
+ that is part of the Library, the object code for the work may be a
+ derivative work of the Library even though the source code is not.
+ Whether this is true is especially significant if the work can be
+ linked without the Library, or if the work is itself a library. The
+ threshold for this to be true is not precisely defined by law.
+
+ If such an object file uses only numerical parameters, data
+ structure layouts and accessors, and small macros and small inline
+ functions (ten lines or less in length), then the use of the object
+ file is unrestricted, regardless of whether it is legally a derivative
+ work. (Executables containing this object code plus portions of the
+ Library will still fall under Section 6.)
+
+ Otherwise, if the work is a derivative of the Library, you may
+ distribute the object code for the work under the terms of Section 6.
+ Any executables containing that work also fall under Section 6,
+ whether or not they are linked directly with the Library itself.
+
+ 6. As an exception to the Sections above, you may also combine or
+ link a "work that uses the Library" with the Library to produce a
+ work containing portions of the Library, and distribute that work
+ under terms of your choice, provided that the terms permit
+ modification of the work for the customer's own use and reverse
+ engineering for debugging such modifications.
+
+ You must give prominent notice with each copy of the work that the
+ Library is used in it and that the Library and its use are covered by
+ this License. You must supply a copy of this License. If the work
+ during execution displays copyright notices, you must include the
+ copyright notice for the Library among them, as well as a reference
+ directing the user to the copy of this License. Also, you must do one
+ of these things:
+
+ a) Accompany the work with the complete corresponding
+ machine-readable source code for the Library including whatever
+ changes were used in the work (which must be distributed under
+ Sections 1 and 2 above); and, if the work is an executable linked
+ with the Library, with the complete machine-readable "work that
+ uses the Library", as object code and/or source code, so that the
+ user can modify the Library and then relink to produce a modified
+ executable containing the modified Library. (It is understood
+ that the user who changes the contents of definitions files in the
+ Library will not necessarily be able to recompile the application
+ to use the modified definitions.)
+
+ b) Use a suitable shared library mechanism for linking with the
+ Library. A suitable mechanism is one that (1) uses at run time a
+ copy of the library already present on the user's computer system,
+ rather than copying library functions into the executable, and (2)
+ will operate properly with a modified version of the library, if
+ the user installs one, as long as the modified version is
+ interface-compatible with the version that the work was made with.
+
+ c) Accompany the work with a written offer, valid for at
+ least three years, to give the same user the materials
+ specified in Subsection 6a, above, for a charge no more
+ than the cost of performing this distribution.
+
+ d) If distribution of the work is made by offering access to copy
+ from a designated place, offer equivalent access to copy the above
+ specified materials from the same place.
+
+ e) Verify that the user has already received a copy of these
+ materials or that you have already sent this user a copy.
+
+ For an executable, the required form of the "work that uses the
+ Library" must include any data and utility programs needed for
+ reproducing the executable from it. However, as a special exception,
+ the materials to be distributed need not include anything that is
+ normally distributed (in either source or binary form) with the major
+ components (compiler, kernel, and so on) of the operating system on
+ which the executable runs, unless that component itself accompanies
+ the executable.
+
+ It may happen that this requirement contradicts the license
+ restrictions of other proprietary libraries that do not normally
+ accompany the operating system. Such a contradiction means you cannot
+ use both them and the Library together in an executable that you
+ distribute.
+
+ 7. You may place library facilities that are a work based on the
+ Library side-by-side in a single library together with other library
+ facilities not covered by this License, and distribute such a combined
+ library, provided that the separate distribution of the work based on
+ the Library and of the other library facilities is otherwise
+ permitted, and provided that you do these two things:
+
+ a) Accompany the combined library with a copy of the same work
+ based on the Library, uncombined with any other library
+ facilities. This must be distributed under the terms of the
+ Sections above.
+
+ b) Give prominent notice with the combined library of the fact
+ that part of it is a work based on the Library, and explaining
+ where to find the accompanying uncombined form of the same work.
+
+ 8. You may not copy, modify, sublicense, link with, or distribute
+ the Library except as expressly provided under this License. Any
+ attempt otherwise to copy, modify, sublicense, link with, or
+ distribute the Library is void, and will automatically terminate your
+ rights under this License. However, parties who have received copies,
+ or rights, from you under this License will not have their licenses
+ terminated so long as such parties remain in full compliance.
+
+ 9. You are not required to accept this License, since you have not
+ signed it. However, nothing else grants you permission to modify or
+ distribute the Library or its derivative works. These actions are
+ prohibited by law if you do not accept this License. Therefore, by
+ modifying or distributing the Library (or any work based on the
+ Library), you indicate your acceptance of this License to do so, and
+ all its terms and conditions for copying, distributing or modifying
+ the Library or works based on it.
+
+ 10. Each time you redistribute the Library (or any work based on the
+ Library), the recipient automatically receives a license from the
+ original licensor to copy, distribute, link with or modify the Library
+ subject to these terms and conditions. You may not impose any further
+ restrictions on the recipients' exercise of the rights granted herein.
+ You are not responsible for enforcing compliance by third parties with
+ this License.
+
+ 11. If, as a consequence of a court judgment or allegation of patent
+ infringement or for any other reason (not limited to patent issues),
+ conditions are imposed on you (whether by court order, agreement or
+ otherwise) that contradict the conditions of this License, they do not
+ excuse you from the conditions of this License. If you cannot
+ distribute so as to satisfy simultaneously your obligations under this
+ License and any other pertinent obligations, then as a consequence you
+ may not distribute the Library at all. For example, if a patent
+ license would not permit royalty-free redistribution of the Library by
+ all those who receive copies directly or indirectly through you, then
+ the only way you could satisfy both it and this License would be to
+ refrain entirely from distribution of the Library.
+
+ If any portion of this section is held invalid or unenforceable under any
+ particular circumstance, the balance of the section is intended to apply,
+ and the section as a whole is intended to apply in other circumstances.
+
+ It is not the purpose of this section to induce you to infringe any
+ patents or other property right claims or to contest validity of any
+ such claims; this section has the sole purpose of protecting the
+ integrity of the free software distribution system which is
+ implemented by public license practices. Many people have made
+ generous contributions to the wide range of software distributed
+ through that system in reliance on consistent application of that
+ system; it is up to the author/donor to decide if he or she is willing
+ to distribute software through any other system and a licensee cannot
+ impose that choice.
+
+ This section is intended to make thoroughly clear what is believed to
+ be a consequence of the rest of this License.
+
+ 12. If the distribution and/or use of the Library is restricted in
+ certain countries either by patents or by copyrighted interfaces, the
+ original copyright holder who places the Library under this License may add
+ an explicit geographical distribution limitation excluding those countries,
+ so that distribution is permitted only in or among countries not thus
+ excluded. In such case, this License incorporates the limitation as if
+ written in the body of this License.
+
+ 13. The Free Software Foundation may publish revised and/or new
+ versions of the Lesser General Public License from time to time.
+ Such new versions will be similar in spirit to the present version,
+ but may differ in detail to address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the Library
+ specifies a version number of this License which applies to it and
+ "any later version", you have the option of following the terms and
+ conditions either of that version or of any later version published by
+ the Free Software Foundation. If the Library does not specify a
+ license version number, you may choose any version ever published by
+ the Free Software Foundation.
+
+ 14. If you wish to incorporate parts of the Library into other free
+ programs whose distribution conditions are incompatible with these,
+ write to the author to ask for permission. For software which is
+ copyrighted by the Free Software Foundation, write to the Free
+ Software Foundation; we sometimes make exceptions for this. Our
+ decision will be guided by the two goals of preserving the free status
+ of all derivatives of our free software and of promoting the sharing
+ and reuse of software generally.
+
+ NO WARRANTY
+
+ 15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO
+ WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW.
+ EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR
+ OTHER PARTIES PROVIDE THE LIBRARY "AS IS" WITHOUT WARRANTY OF ANY
+ KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE
+ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE
+ LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME
+ THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN
+ WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY
+ AND/OR REDISTRIBUTE THE LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU
+ FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR
+ CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE
+ LIBRARY (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING
+ RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A
+ FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF
+ SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH
+ DAMAGES.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Libraries
+
+ If you develop a new library, and you want it to be of the greatest
+ possible use to the public, we recommend making it free software that
+ everyone can redistribute and change. You can do so by permitting
+ redistribution under these terms (or, alternatively, under the terms of the
+ ordinary General Public License).
+
+ To apply these terms, attach the following notices to the library. It is
+ safest to attach them to the start of each source file to most effectively
+ convey the exclusion of warranty; and each file should have at least the
+ "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the library's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ This library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with this library; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+
+ Also add information on how to contact you by electronic and paper mail.
+
+ You should also get your employer (if you work as a programmer) or your
+ school, if any, to sign a "copyright disclaimer" for the library, if
+ necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the
+ library `Frob' (a library for tweaking knobs) written by James Random Hacker.
+
+ <signature of Ty Coon>, 1 April 1990
+ Ty Coon, President of Vice
+
+ That's all there is to it!
+
All deal.II users (and everybody interested)
are encouraged to join the deal.II
- Discussion group on Google Groups.
+ Discussion Group on Google Groups.
@@ -44,7 +44,7 @@
the
the deal.II
Discussion Group page and click on the corresponding
- button (preferred), or you can send email
+ button, or you can send email
to dealii@googlegroups.com. Or
just post from the window below. Note that you have to be a
group member with your specific sender address in order to be
@@ -71,16 +71,16 @@
flaming them on missing items) should be sent to one of
- authors at dealii.org
+ TODO
- Wolfgang.Bangerth at dealii.org
+ Wolfgang Bangerth < bangerth@gmail.com >
- Timo.Heister at dealii.org
+ Timo Heister < heister@math.tamu.edu >
- Guido.Kanschat at dealii.org
+ Guido Kanschat < dr.guido.kanschat@gmail.com >
Everything, that is not directly related to the library is to be considered
diff --git a/deal.II/doc/navbar.html b/deal.II/doc/navbar.html
index 64b987fe61..fcc7b54ec1 100644
--- a/deal.II/doc/navbar.html
+++ b/deal.II/doc/navbar.html
@@ -13,13 +13,22 @@
Development Publications
diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 399aac68cd..a29684ee60 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -64,6 +64,26 @@ this function.
+
+
New: All vector classes now have a member function
+ locally_owned_elements that returns an index
+ set indicating which elements of this vector the current
+ processor owns.
+
+ (Wolfgang Bangerth, 2013/05/24)
+
+
+
+
New: A new element FE_Q_iso_Q1 has been implemented that is defined by
+ a subdivision of the element into smaller Q1 elements. An element of order
+ @p p is similar to FE_Q of degree @p p with the same numbering of degrees of
+ freedom. The element is useful e.g. for defining a sparser preconditioner
+ matrix for AMG at higher order FE_Q elements or for representing a component
+ of a system of PDEs where higher resolution is preferred over high order.
+
+ (Martin Kronbichler, 2013/05/14)
+
+
New: The step-49 tutorial program now also has a discussion on
what to do once you have a coarse mesh and want to refine it.
@@ -110,6 +130,152 @@ this function.
Specific improvements
+
New: ConstraintMatrix::distribute is now also implemented for
+arguments of type PETScWrappers::MPI::BlockVector.
+
+(Wolfgang Bangerth, 2013/05/25)
+
+
+
Fixed: IndexSet::operator== returned the wrong results
+in some cases.
+
+(Wolfgang Bangerth, 2013/05/25)
+
+
+
New: The global function complete_index_set()
+creates and returns an index set of given size that contains
+every single index with this range.
+
+(Wolfgang Bangerth, 2013/05/24)
+
+
+
New: All vector classes now have a static member variable
+supports_distributed_data that indicates whether the
+vector class supports data that is distributed across multiple
+processors. This variable is provided as a traits variable
+to allow generic algorithms working on general vector types to
+query the capabilities of the vector class at compile time.
+
+(Wolfgang Bangerth, 2013/05/23)
+
+
+
Fixed: FETools::back_interpolate has been revised to work correctly
+also with parallel::distributed::Vector.
+
+(Martin Steigemann, 2013/05/23)
+
+
+
Removed: The file mesh_worker/vector_info.h was unused and
+untested. It has thus been removed.
+
+(Wolfgang Bangerth, Guido Kanschat, 2013/05/21)
+
+
+
Fixed: The method parallel::distributed::Vector::compress
+(VectorOperation::insert) previously set the elements of ghost elements
+unconditionally on the owning processor, even if they had not been touched.
+This led to a problem in certain library functions where vector entries became
+zero in a spurious way. This is now fixed by discarding the elements in ghost
+entries for the VectorOperation::insert operation. This is legitimate since we
+assume consistency of set elements across processors, so the owning processor
+sets the element already.
+
+(Martin Kronbichler, 2013/05/21)
+
+
+
Improved: DoFTools::make_periodicity_constraints now also works
+for meshes where the refinement level of the two sides of the domain
+is not the same, i.e., one side is more refined than the other.
+
+(Wolfgang Bangerth, 2013/05/20)
+
+
+
Improved: Through the fields DataOutBase::VtkFlags::time and
+DataOutBase::VtkFlags::cycle, it is now possible to encode the time and/or
+cycle within a nonlinear or other iteration in VTK and VTU files written
+via DataOutBase::write_vtk and DataOutBase::write_vtu.
+
+(Wolfgang Bangerth, 2013/05/12)
+
+
+
Fixed: The method ConvergenceTable::evaluate_convergence_rates with
+ reference column did not take the dimension of the reference column into
+ account, leading to wrong logarithmic rates for dim!=2. This can now be fixed
+ by specifying the dimension as a last argument.
+
+(Martin Kronbichler, 2013/05/10)
+
+
+
Improved: The functions MatrixTools::create_mass_matrix and
+MatrixTools::create_laplace_matrix take now an optional ConstraintMatrix
+argument that allows to directly apply the constraints. This also helps
+VectorTools::project. Note that not providing constraints remains the default
+and recommended way to ensure consistency when several matrices are combined.
+
+(Martin Kronbichler, 2013/05/08)
+
+
+
New: The classes TrilinosWrappers::SparseMatrix and
+TrilinosWrappers::BlockSparseMatrix now fully implement vmult and Tvmult with
+deal.II's own vector classes Vector and
+parallel::distributed::Vector.
+
+(Martin Kronbichler, 2013/05/08)
+
+
+
Improved: The matrix-vector product ChunkSparseMatrix::vmult now runs in
+parallel in shared memory.
+
+(Martin Kronbichler, 2013/05/07)
+
+
+
New: The class ChunkSparseMatrix and the associated
+ChunkSparsityPattern now offer iterator classes to iterate over rows of the
+matrix in an STL-like way.
+
+(Martin Kronbichler, 2013/05/07)
+
+
+
Fixed: The stopping criterion for early exit in SolverBicgstab did not
+work properly for systems with large values, leading to premature exit. This
+is now fixed.
+
+(Martin Kronbichler, 2013/05/07)
+
+
+
Changed: The SolverGMRES implementation previously applied two
+iterations of the modified Gram–Schmidt algorithm for
+orthogonalization. In many situations one iteration is enough. The algorithm
+can now detect loss of orthogonality and enables re-orthogonalization only if
+necessary. The second iteration (and, hence, old behavior) can be forced by
+the flag SolverGMRES::AdditionalData::force_re_orthogonalization.
+
+(Martin Kronbichler, 2013/05/06)
+
+
+
Changed: FETools::interpolate is now instantiated for all
+vector types, not just dealii::Vector and dealii::BlockVector.
+
+(Wolfgang Bangerth, 2013/05/06)
+
+
+
Fixed: setting values in TrilinosWrappers::SparseMatrix
+in parallel was adding the values instead.
+
+(Bruno Turcksin, Timo Heister, 2013/05/03)
+
+
+
Fixed: Generate an error if the user tries to refine a cell
+that is already on the maximum level in a distributed triangulation.
+
+(Timo Heister, 2013/05/01)
+
+
+
Fixed: The version of ParameterHandler::set that takes a boolean
+as second argument was broken and did not work. This is now fixed.
+
+(Ashkan Dorostkar, Wolfgang Bangerth, 2013/04/30)
+
Fixed: PETScWrappers::VectorBase::print now saves and restores
the precision
@@ -261,13 +427,17 @@ sense.
(Guido Kanschat, 2013/03/21)
-
Added GridOut::write_svg to allow for the output of two-dimensional
-triangulations in two space dimensions in the SVG format (Scalable Vector
-Graphics, an XML-based vector image format recommended by the World
-Wide Web Consortium W3C). This function also provides cell coloring
-and cell labeling for the visualization of basic cell properties.
+
Added GridOut::write_svg() to allow for the output of
+two-dimensional triangulations in two space dimensions in the SVG
+format (Scalable Vector Graphics, an generic XML-based vector image
+format developed and maintained by the World Wide Web Consortium W3C).
+This function also provides cell coloring and cell labeling for the
+visualization of basic cell properties. Pespective view is further
+possible and the cell level number may be converted into altitude,
+revealing the inactive cells lying below.
(Christian Wülker, 2013/03/21)
+
Added TimerOutput::reset to remove the collected information so far and
added a new frequency TimerOutput::never to only output information if
@@ -292,6 +462,16 @@ This is now fixed.
(Timo Heister, 2013/03/01)
+
Added DataOutBase::write_svg() to allow for the output of a given
+list of patches in two space dimensions in the SVG format (Scalable Vector
+Graphics, an generic XML-based vector image format developed and maintained
+by the World Wide Web Consortium W3C). An additional dimension (z-direction)
+is employed for the visualization of data values taken from a data vector.
+This function also provides patch coloring for the visual enhancement.
+
+(Christian Wülker, 2013/05/10)
+
+
diff --git a/deal.II/doc/publications/index.html b/deal.II/doc/publications/index.html
index 346587c363..77d0a3a83b 100644
--- a/deal.II/doc/publications/index.html
+++ b/deal.II/doc/publications/index.html
@@ -124,7 +124,7 @@
@Manual{DealIIReference,
title = {{\tt deal.{I}{I}} Differential Equations Analysis Library,
Technical Reference},
- author = {W. Bangerth and T. Heister and G. Kanschat},
+ author = {W. Bangerth and T. Heister and G. Kanschat and others},
note = {\texttt{http://www.dealii.org}},
url = {http://www.dealii.org}
}
@@ -297,6 +297,15 @@
+
K. Dugan
+
+ Dynamic adaptive multimesh refinement for coupled
+ physics applicable to nuclear engineering
+
+
+ MSc thesis, Texas A&M University, 2013.
+
+
Y. Efendiev, J. Galvis, S. Pauletti
Multiscale finite element methods for flows on rough
@@ -373,6 +382,13 @@
Series, pp. 79-86, 2013.
+
J. Reinhardt, F. Weysser, J. M. Brader
+
+ Density functional approach to nonlinear rheology
+
+ Europhysics Letters, vol. 102, article 28011, 2013.
+
+
S. Rudraraju, K. L. Mills, R. Kemkemer, K. Garikipati
Multiphysics Modeling of Reactions, Mass Transport and Mechanics of Tumor Growth
@@ -382,6 +398,15 @@
(eds.), Spinger Verlag, pp. 293-303, 2013.
+
N. Vakulin, R. Shaw, T. Livingston
+
+ ELEC 490, final report: High performance computing with GPUs
+
+
+ Semester project for ELEC 490, final report. Queen's University,
+ 2013
+
+
S. Yoon, Y. T. Kang
Mass transfer enhancement in non-Brownian particle
diff --git a/deal.II/doc/readme.html b/deal.II/doc/readme.html
index 6663f289e0..e80ebdac83 100644
--- a/deal.II/doc/readme.html
+++ b/deal.II/doc/readme.html
@@ -785,50 +785,59 @@
License
- We have placed deal.II under an Open Source
- license, which allows you to use the library free of charge. You
- are guaranteed full access to the source code and are encouraged
- to help in the further development of the library. Follow this link to read the full
- text of the license,
+ The deal.II library is free software; you can use it, redistribute
+ it, and/or modify it under the terms of the
+ GNU Lesser General
+ Public License (LGPL) as published by the Free Software
+ Foundation; either version 2.1 of the License, or (at your option)
+ any later version.
+
+
+ This allows you to use the library free of charge for private,
+ academic, or commercial use (under the terms of the LGPL v2.1 or
+ later). You are guaranteed full access to the source code and are
+ encouraged to help in the further development of the library. Follow
+ this link for detailed
+ information.
-
- The basics of the license are in short:
+ Please note:
- The library comes free of charge for all users, but
- we keep a copyright in it. We will grant everyone access to
- the library and its source code, provided they keep to the
- license.
-
+ Detailed license information can be found following
+ this link.
- The libraries may not be stripped of our copyright. We would
- appreciate being contacted if you want to use the library
- commercially.
-
+ As a contributor to this project, you agree that any of your
+ contributions be licensed under the same terms and conditions as
+ the license of the deal.II project granted to you.
+
+ We, the deal.II authors, do not require
+ copyright assignments for contributions. This means that the
+ copyright for code contributions in the deal.II project is held by
+ its respective contributors who have each agreed to release their
+ contributed code under the terms of the LGPL v2.1 or later.
+
+ In addition to the terms imposed by the LGPL (version 2.1 or
+ later), we ask for the courtesy that scientific publications
+ presenting results obtained with this libraries acknowledge its
+ use. Please cite one of the papers referenced
+ here.
- We ask that scientific publications presenting results obtained with this
- libraries acknowledge its use. Please cite one of the papers referenced
- here.
+ deal.II can interface with a number of other packages that you either
+ have to install yourself. They are, of course, covered by their
+ own licenses. In addition, deal.II comes bundled with copies of
+ UMFPACK,
+ Threading Building Blocks,
+ BOOST and
+ FunctionParser, courtesy of their authors.
+ These are also covered by their own licenses; please refer to
+ their webpages for more information.
-
- deal.II can interface with a number of
- other packages that you have to install
- yourself. They are, of course, covered by their own licenses. In addition,
- deal.II comes bundled with copies of
- UMFPACK,
- Threading Building Blocks,
- BOOST
- and
- FunctionParser,
- courtesy of their authors. These are also covered by their
- own licenses; please refer to their webpages for more information.
-
diff --git a/deal.II/examples/step-26/step-26.cc b/deal.II/examples/step-26/step-26.cc
index a4528ae963..0d4da01968 100644
--- a/deal.II/examples/step-26/step-26.cc
+++ b/deal.II/examples/step-26/step-26.cc
@@ -15,6 +15,7 @@
#include
#include
#include
+#include
#include
#include
#include
@@ -172,8 +173,6 @@ namespace Step26
template
void HeatEquation::setup_system()
{
-//TODO: we need to take care of constraints here! :-)
-
dof_handler.distribute_dofs(fe);
std::cout << std::endl
@@ -185,12 +184,17 @@ namespace Step26
<< std::endl
<< std::endl;
- sparsity_pattern.reinit(dof_handler.n_dofs(),
- dof_handler.n_dofs(),
- dof_handler.max_couplings_between_dofs());
+ constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+ constraints.close();
+
+ CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
- sparsity_pattern);
- sparsity_pattern.compress();
+ c_sparsity,
+ constraints,
+ /*keep_constrained_dofs = */ true);
+ sparsity_pattern.copy_from(c_sparsity);
mass_matrix.reinit(sparsity_pattern);
laplace_matrix.reinit(sparsity_pattern);
@@ -198,16 +202,18 @@ namespace Step26
MatrixCreator::create_mass_matrix(dof_handler,
QGauss(fe.degree+1),
- mass_matrix);
+ mass_matrix,
+ (const Function*)0,
+ constraints);
MatrixCreator::create_laplace_matrix(dof_handler,
QGauss(fe.degree+1),
- laplace_matrix);
+ laplace_matrix,
+ (const Function*)0,
+ constraints);
solution.reinit(dof_handler.n_dofs());
old_solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
-
- constraints.close();
}
@@ -223,6 +229,8 @@ namespace Step26
cg.solve(system_matrix, solution, system_rhs, preconditioner);
+ constraints.distribute(solution);
+
std::cout << " " << solver_control.last_step()
<< " CG iterations." << std::endl;
}
@@ -362,8 +370,8 @@ namespace Step26
template
void HeatEquation::run()
{
- const unsigned int initial_global_refinement = (dim == 2 ? 4 : 2);
- const unsigned int n_adaptive_pre_refinement_steps = 3;
+ const unsigned int initial_global_refinement = (dim == 2 ? 1 : 2);
+ const unsigned int n_adaptive_pre_refinement_steps = 1;
GridGenerator::hyper_L (triangulation);
triangulation.refine_global (initial_global_refinement);
@@ -440,6 +448,8 @@ start_time_iteration:
system_rhs);
}
+ constraints.condense (system_rhs);
+
solve_time_step();
output_results();
diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc
index de28338a59..54de0ab8f4 100644
--- a/deal.II/examples/step-42/step-42.cc
+++ b/deal.II/examples/step-42/step-42.cc
@@ -1551,7 +1551,7 @@ namespace Step42
{
pcout << "Read the obstacle from a file." << std::endl;
input_obstacle.reset (new Input("obstacle_file.pbm"));
- pcout << "Ostacle is available now." << std::endl;
+ pcout << "Obstacle is available now." << std::endl;
}
const unsigned int n_cycles = 6;
diff --git a/deal.II/examples/step-48/step-48.cc b/deal.II/examples/step-48/step-48.cc
index e5e71e7d20..79250490ed 100644
--- a/deal.II/examples/step-48/step-48.cc
+++ b/deal.II/examples/step-48/step-48.cc
@@ -3,7 +3,7 @@
/* $Id$ */
/* */
-/* Copyright (C) 2011-2012 by the deal.II authors */
+/* Copyright (C) 2011-2013 by the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
@@ -603,13 +603,20 @@ namespace Step48
// @sect3{The main function}
-// This is as in all other programs:
+// As in step-40, we initialize MPI at the start of the program. Since we will
+// in general mix MPI parallelization with threads, we also set the third
+// argument in MPI_InitFinalize that controls the number of threads to an
+// invalid number, which means that the TBB library chooses the number of
+// threads automatically, typically to the number of available cores in the
+// system. As an alternative, you can also set this number manually if you
+// want to set a specific number of threads (e.g. when MPI-only is required).
int main (int argc, char **argv)
{
using namespace Step48;
using namespace dealii;
- Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
+ Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv,
+ numbers::invalid_unsigned_int);
try
{
diff --git a/deal.II/examples/step-50/CMakeLists.txt b/deal.II/examples/step-50/CMakeLists.txt
new file mode 100644
index 0000000000..babcd817e0
--- /dev/null
+++ b/deal.II/examples/step-50/CMakeLists.txt
@@ -0,0 +1,52 @@
+##
+# CMake script for the step-32 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-50")
+
+# Declare all source files the target consists of:
+SET(TARGET_SRC
+ ${TARGET}.cc
+ # You can specify additional files here!
+ )
+
+# A custom command line to run the program
+SET(TARGET_RUN mpirun -np 2 ${TARGET})
+
+# Define the output that should be cleaned:
+SET(CLEAN_UP_FILES *.vtu *.pvtu *.visit)
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
+
+FIND_PACKAGE(deal.II 8.0 QUIET
+ HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+ )
+IF(NOT ${deal.II_FOUND})
+ MESSAGE(FATAL_ERROR "\n"
+ "*** Could not locate deal.II. ***\n\n"
+ "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+ "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+ )
+ENDIF()
+
+#
+# Are all dependencies fullfilled?
+#
+IF( NOT DEAL_II_WITH_MPI OR
+ NOT DEAL_II_WITH_P4EST OR
+ NOT DEAL_II_WITH_TRILINOS )
+ MESSAGE(FATAL_ERROR "
+Error! The deal.II library found at ${DEAL_II_PATH} was not configured with
+ DEAL_II_WITH_MPI = ON
+ DEAL_II_WITH_P4EST = ON
+ DEAL_II_WITH_TRILINOS = ON
+which is required for this tutorial step."
+ )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/deal.II/examples/step-50/step-50.cc b/deal.II/examples/step-50/step-50.cc
index 640c4b38be..abca04f105 100644
--- a/deal.II/examples/step-50/step-50.cc
+++ b/deal.II/examples/step-50/step-50.cc
@@ -11,6 +11,8 @@
/* to the file deal.II/doc/license.html for the text and */
/* further information on this license. */
+// parallel geometric multi-grid. work in progress!
+
// As discussed in the introduction, most of
// this program is copied almost verbatim
// from step-6, which itself is only a slight
@@ -47,6 +49,7 @@
#include
#include
#include
+#include
#include
#include
@@ -132,7 +135,6 @@ namespace Step50
typedef TrilinosWrappers::SparseMatrix matrix_t;
typedef TrilinosWrappers::MPI::Vector vector_t;
- SparsityPattern sparsity_pattern;
matrix_t system_matrix;
// We need an additional object for the
@@ -188,7 +190,6 @@ namespace Step50
// transpose of the other, and so we only
// have to build one; we choose the one
// from coarse to fine.
- MGLevelObject mg_sparsity_patterns;
MGLevelObject mg_matrices;
MGLevelObject mg_interface_matrices;
MGConstrainedDoFs mg_constrained_dofs;
@@ -279,7 +280,8 @@ namespace Step50
LaplaceProblem::LaplaceProblem (const unsigned int degree)
:
triangulation (MPI_COMM_WORLD,Triangulation::
- limit_level_difference_at_vertices),
+ limit_level_difference_at_vertices,
+ parallel::distributed::Triangulation::construct_multigrid_hierarchy),
fe (degree),
mg_dof_handler (triangulation),
degree(degree)
@@ -296,6 +298,21 @@ namespace Step50
template
void LaplaceProblem::setup_system ()
{
+ std::string filename = "grid-"
+ + Utilities::int_to_string(triangulation.n_levels(), 2)
+ + "-"
+ + Utilities::int_to_string (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD), 3)
+ + ".svg";
+
+ std::ofstream out (filename.c_str());
+ GridOut grid_out;
+ GridOutFlags::Svg svg_flags;
+ svg_flags.polar_angle = 60;
+ svg_flags.label_level_subdomain_id = true;
+ svg_flags.convert_level_number_to_height = true;
+ grid_out.set_flags(svg_flags);
+ grid_out.write_svg (triangulation, out);
+
mg_dof_handler.distribute_dofs (fe);
mg_dof_handler.distribute_mg_dofs (fe);
@@ -312,11 +329,6 @@ namespace Step50
<< mg_dof_handler.n_dofs(l);
deallog << std::endl;
- sparsity_pattern.reinit (mg_dof_handler.n_dofs(),
- mg_dof_handler.n_dofs(),
- mg_dof_handler.max_couplings_between_dofs());
- DoFTools::make_sparsity_pattern (mg_dof_handler, sparsity_pattern);
-
//solution.reinit (mg_dof_handler.n_dofs());
//system_rhs.reinit (mg_dof_handler.n_dofs());
@@ -356,9 +368,10 @@ namespace Step50
constraints);
constraints.close ();
hanging_node_constraints.close ();
- constraints.condense (sparsity_pattern);
- sparsity_pattern.compress();
- system_matrix.reinit (sparsity_pattern);
+
+ CompressedSimpleSparsityPattern csp(mg_dof_handler.n_dofs(), mg_dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern (mg_dof_handler, csp, constraints);
+ system_matrix.reinit (mg_dof_handler.locally_owned_dofs(), csp, MPI_COMM_WORLD, true);
// The multigrid constraints have to be
// initialized. They need to know about
@@ -390,7 +403,6 @@ namespace Step50
mg_interface_matrices.clear ();
mg_matrices.resize(0, n_levels-1);
mg_matrices.clear ();
- mg_sparsity_patterns.resize(0, n_levels-1);
// Now, we have to provide a matrix on each
// level. To this end, we first use the
@@ -425,10 +437,10 @@ namespace Step50
mg_dof_handler.n_dofs(level));
MGTools::make_sparsity_pattern(mg_dof_handler, csp, level);
- mg_sparsity_patterns[level].copy_from (csp);
-
- mg_matrices[level].reinit(mg_sparsity_patterns[level]);
- mg_interface_matrices[level].reinit(mg_sparsity_patterns[level]);
+ mg_matrices[level].reinit(mg_dof_handler.locally_owned_mg_dofs(level),
+ mg_dof_handler.locally_owned_mg_dofs(level),
+ csp,
+ MPI_COMM_WORLD, true);
}
}
@@ -971,12 +983,12 @@ namespace Step50
if (cycle == 0)
{
- GridGenerator::hyper_ball (triangulation);
+ GridGenerator::hyper_cube (triangulation);
- static const HyperBallBoundary boundary;
- triangulation.set_boundary (0, boundary);
+ // static const HyperBallBoundary boundary;
+ // triangulation.set_boundary (0, boundary);
- triangulation.refine_global (1);
+ triangulation.refine_global (2);
}
else
refine_grid ();
diff --git a/deal.II/include/deal.II/base/convergence_table.h b/deal.II/include/deal.II/base/convergence_table.h
index df3bed1a50..a262a6be3b 100644
--- a/deal.II/include/deal.II/base/convergence_table.h
+++ b/deal.II/include/deal.II/base/convergence_table.h
@@ -73,175 +73,127 @@ public:
*/
none,
/**
- * Quotient of values in
- * the previous row and in
- * this row.
+ * Quotient of values in the previous row and in this row.
*/
reduction_rate,
/**
- * Logarithm of
- * #reduction_rate to the
- * base 2 representing the
- * order of convergence
- * when halving the grid
- * size, e.g. from h to
- * h/2.
+ * Logarithm of #reduction_rate to the base 2 representing the order of
+ * convergence when halving the grid size, e.g. from h to h/2.
*/
reduction_rate_log2
};
/**
- * Evaluates the convergence rates of the
- * data column data_column_key
- * due to the #RateMode in relation to
- * the reference column reference_column_key.
- * Be sure that the value types of the
- * table entries of the
- * data column and the reference data column
- * is a number, i.e. double, float,
- * (unsigned) int, and so on.
+ * Evaluates the convergence rates of the data column
+ * data_column_key due to the #RateMode in relation to the
+ * reference column reference_column_key. Be sure that the value
+ * types of the table entries of the data column and the reference data
+ * column is a number, i.e. double, float, (unsigned) int, and so on.
*
- * The new rate column and the data column
- * will be merged to a supercolumn. The
- * tex caption of the supercolumn will be
- * (by default) the same as the one of the
- * data column. This may be changed by using
- * the set_tex_supercaption (...) function
- * of the base class TableHandler.
+ * As this class has no information on the space dimension upon which the
+ * reference column vs. the value column is based upon, it needs to be
+ * passed as last argument to this method. The default dimension for the
+ * reference column is 2, which is appropriate for the number of cells
+ * in 2D. If you work in 3D, set the number to 3. If the reference column is
+ * $1/h$, remember to set the dimension to 1 also when working in 3D to get
+ * correct rates.
*
- * This method behaves in the following way:
+ * The new rate column and the data column will be merged to a
+ * supercolumn. The tex caption of the supercolumn will be (by default) the
+ * same as the one of the data column. This may be changed by using the
+ * set_tex_supercaption (...) function of the base class
+ * TableHandler.
*
- * If RateMode is reduction_rate, then the computed
- * output is
- * $ \frac{e_{n-1}/k_{n-1}}{e_n/k_n} $.
+ * This method behaves in the following way:
*
- * Where $k$ is the reference column.
+ * If RateMode is reduction_rate, then the computed output is
+ * $ \frac{e_{n-1}/k_{n-1}}{e_n/k_n}, $
+ * where $k$ is the reference column (no dimension dependence!).
*
- * If RateMode is reduction_rate_log2, then the
- * computed output is
+ * If RateMode is reduction_rate_log2, then the computed output is
* $
- * 2\frac{\log |e_{n-1}/e_{n}|}{\log |k_n/k_{n-1}|}
+ * dim \frac{\log |e_{n-1}/e_{n}|}{\log |k_n/k_{n-1}|}
* $.
*
- * This is useful, for example, if we use as
- * reference key the number of degrees of freedom.
- * Assuming that the error is proportional to
- * $ C (1/\sqrt{k})^r $, then this method will
- * produce the rate $r$ as a result.
+ * This is useful, for example, if we use as reference key the number of
+ * degrees of freedom or better, the number of cells. Assuming that the
+ * error is proportional to $ C (1/\sqrt{k})^r $ in 2D, then this method
+ * will produce the rate $r$ as a result. For general dimension, as described
+ * by the last parameter of this function, the formula needs to be
+ * $ C (1/\sqrt[dim]{k})^r $.
*
- * @note Since this function adds columns
- * to the table after several rows have
- * already been filled, it switches off
- * the auto fill mode of the TableHandler
- * base class. If you intend to add
- * further data with auto fill, you will
- * have to re-enable it after calling
- * this function.
+ * @note Since this function adds columns to the table after several rows
+ * have already been filled, it switches off the auto fill mode of the
+ * TableHandler base class. If you intend to add further data with auto
+ * fill, you will have to re-enable it after calling this function.
*/
void
evaluate_convergence_rates (const std::string &data_column_key,
const std::string &reference_column_key,
- const RateMode rate_mode);
+ const RateMode rate_mode,
+ const unsigned int dim = 2);
/**
- * Evaluates the convergence rates of the
- * data column data_column_key
- * due to the #RateMode.
- * Be sure that the value types of the
- * table entries of the data column
- * is a number, i.e. double, float,
- * (unsigned) int, and so on.
+ * Evaluates the convergence rates of the data column
+ * data_column_key due to the #RateMode. Be sure that the value
+ * types of the table entries of the data column is a number, i.e. double,
+ * float, (unsigned) int, and so on.
*
- * The new rate column and the data column
- * will be merged to a supercolumn. The
- * tex caption of the supercolumn will be
- * (by default) the same as the one of the
- * data column. This may be changed by using
- * the set_tex_supercaption() function
- * of the base class TableHandler.
+ * The new rate column and the data column will be merged to a
+ * supercolumn. The tex caption of the supercolumn will be (by default) the
+ * same as the one of the data column. This may be changed by using the
+ * set_tex_supercaption() function of the base class TableHandler.
*
- * @note Since this function adds columns
- * to the table after several rows have
- * already been filled, it switches off
- * the auto fill mode of the TableHandler
- * base class. If you intend to add
- * further data with auto fill, you will
- * have to re-enable it after calling
- * this function.
+ * @note Since this function adds columns to the table after several rows
+ * have already been filled, it switches off the auto fill mode of the
+ * TableHandler base class. If you intend to add further data with auto
+ * fill, you will have to re-enable it after calling this function.
*/
void
evaluate_convergence_rates (const std::string &data_column_key,
const RateMode rate_mode);
/**
- * Omit this column key
- * (not supercolumn!) from the
- * evaluation of the convergence rates
- * of `all' columns (see the following
- * two functions).
+ * Omit this column key (not supercolumn!) from the evaluation of
+ * the convergence rates of `all' columns (see the following two functions).
*
- * The Column::flag==1 is reserved for
- * omitting the column from convergence
+ * The Column::flag==1 is reserved for omitting the column from convergence
* rate evalution.
*/
void
omit_column_from_convergence_rate_evaluation(const std::string &key);
/**
- * Evaluates convergence rates
- * due to the rate_mode
- * in relation to the reference
- * column
- * reference_column_key. This
- * function evaluates the rates
- * of ALL columns except of the
- * columns that are to be omitted
- * (see previous function) and
- * execpt of the columns that are
- * previously evaluated rate
- * columns. This function allows
- * to evaluate the convergence
- * rate for almost all columns of
- * a table without calling
- * evaluate_convergence_rates()
- * for each column separately.
+ * Evaluates convergence rates due to the rate_mode in relation to
+ * the reference column reference_column_key. This function
+ * evaluates the rates of ALL columns except of the columns that are to be
+ * omitted (see previous function) and execpt of the columns that are
+ * previously evaluated rate columns. This function allows to evaluate the
+ * convergence rate for almost all columns of a table without calling
+ * evaluate_convergence_rates() for each column separately.
*
* Example:
- * Columns like n cells or
- * n dofs columns may be wanted
- * to be omitted in the evaluation
- * of the convergence rates. Hence they
- * should omitted by calling the
- * omit_column_from_convergence_rate_evaluation().
+ * Columns like n cells or n dofs columns may be wanted to
+ * be omitted in the evaluation of the convergence rates. Hence they should
+ * omitted by calling the omit_column_from_convergence_rate_evaluation().
*/
void
evaluate_all_convergence_rates(const std::string &reference_column_key,
const RateMode rate_mode);
/**
- * Evaluates convergence rates
- * due to the rate_mode. This
- * function evaluates the rates of
- * ALL columns except of the
- * columns that are to be omitted
- * (see previous function)
- * and execpt of the columns that are
- * previously
- * evaluated rate columns.
- * This function allows to evaluate
- * the convergence rate for almost all
- * columns of a table without calling
- * evaluate_convergence_rates()
- * for each column seperately.
+ * Evaluates convergence rates due to the rate_mode. This function
+ * evaluates the rates of ALL columns except of the columns that are to be
+ * omitted (see previous function) and execpt of the columns that are
+ * previously evaluated rate columns. This function allows to evaluate the
+ * convergence rate for almost all columns of a table without calling
+ * evaluate_convergence_rates() for each column seperately.
*
* Example:
- * Columns like n cells or
- * n dofs columns may be wanted
- * to be omitted in the evaluation
- * of the convergence rates. Hence they
- * should omitted by calling the
- * omit_column_from_convergence_rate_evaluation().
+ * Columns like n cells or n dofs columns may be wanted to
+ * be omitted in the evaluation of the convergence rates. Hence they should
+ * omitted by calling the omit_column_from_convergence_rate_evaluation().
*/
void
evaluate_all_convergence_rates(const RateMode rate_mode);
diff --git a/deal.II/include/deal.II/base/data_out_base.h b/deal.II/include/deal.II/base/data_out_base.h
index a8124528fc..0ce000d902 100644
--- a/deal.II/include/deal.II/base/data_out_base.h
+++ b/deal.II/include/deal.II/base/data_out_base.h
@@ -21,6 +21,7 @@
#include
#include
+#include
#include
@@ -1186,34 +1187,54 @@ public:
/**
* Flags controlling the details
* of output in VTK
- * format. At present no flags
- * are implemented.
+ * format.
*
* @ingroup output
*/
struct VtkFlags
{
- private:
+ public:
/**
- * Dummy entry to suppress compiler
- * warnings when copying an empty
- * structure. Remove this member
- * when adding the first flag to
- * this structure (and remove the
- * private as well).
+ * The time of the time step if this file is part of a
+ * time dependent simulation.
+ *
+ * The value of this variable is written into the output file according
+ * to the instructions provided in
+ * http://www.visitusers.org/index.php?title=Time_and_Cycle_in_VTK_files
+ * unless it is at its default value of std::numeric_limits::min().
*/
- int dummy;
+ double time;
+
+ /**
+ * The number of the time step if this file is part of a
+ * time dependent simulation, or the cycle within a nonlinear or other
+ * iteration.
+ *
+ * The value of this variable is written into the output file according
+ * to the instructions provided in
+ * http://www.visitusers.org/index.php?title=Time_and_Cycle_in_VTK_files
+ * unless it is at its default value of std::numeric_limits::min().
+ */
+ unsigned int cycle;
- public:
/**
* Default constructor.
*/
- VtkFlags ();
+ VtkFlags (const double time = std::numeric_limits::min(),
+ const unsigned int cycle = std::numeric_limits::min());
/**
- * Declare all flags with name
+ * Declare the flags with name
* and type as offered by this
* class, for use in input files.
+ *
+ * Unlike the flags in many of the other classes similar to this one, we do
+ * not actually declare parameters for the #cycle and #time member variables
+ * of this class. The reason is that there wouldn't appear to be a case where
+ * one would want to declare these parameters in an input file. Rather, these
+ * are typically values that change during the course of a simulation and
+ * can only reasonably be set as part of the execution of a program, rather
+ * than a priori by a user who runs this program.
*/
static void declare_parameters (ParameterHandler &prm);
@@ -1246,6 +1267,77 @@ public:
*/
std::size_t memory_consumption () const;
};
+
+
+ /**
+ * Flags for SVG output.
+ */
+ struct SvgFlags
+ {
+ public:
+ /**
+ * This denotes the number of the
+ * data vector which shall be used
+ * for generating the height
+ * information. By default, the
+ * first data vector is taken,
+ * i.e. height_vector==0, if
+ * there is any data vector. If there
+ * is no data vector, no height
+ * information is generated.
+ */
+ unsigned int height_vector;
+
+ /*
+ * Angles for the perspective view
+ */
+ int azimuth_angle, polar_angle;
+
+ unsigned int line_thickness;
+
+ /*
+ * Draw a margin of 5% around the plotted area
+ */
+ bool margin;
+
+ /*
+ * Draw a colorbar encoding the cell coloring
+ */
+ bool draw_colorbar;
+
+ /*
+ * Constructor.
+ */
+ SvgFlags(const unsigned int height_vector = 0,
+ const int azimuth_angle = 37,
+ const int polar_angle = 45,
+ const unsigned int line_thickness = 1,
+ const bool margin = true,
+ const bool draw_colorbar = true);
+
+ /**
+ * Determine an estimate for
+ * the memory consumption (in
+ * bytes) of this
+ * object. Since sometimes
+ * the size of objects can
+ * not be determined exactly
+ * (for example: what is the
+ * memory consumption of an
+ * STL std::map type with a
+ * certain number of
+ * elements?), this is only
+ * an estimate. however often
+ * quite close to the true
+ * value.
+ */
+ std::size_t memory_consumption () const;
+
+ private:
+
+ };
+
+
/**
* Flags controlling the details
* of output in deal.II
@@ -1398,6 +1490,12 @@ public:
*/
vtu,
+ /**
+ * Output in
+ * SVG format.
+ */
+ svg,
+
/**
* Output in deal.II
* intermediate format.
@@ -1693,6 +1791,13 @@ public:
* provisions that allow these components to be output by a single
* name rather than having to group several scalar fields into a
* vector later on in the visualization program.
+ *
+ * @note VTK is a legacy format and has largely been supplanted by the VTU
+ * format (an XML-structured version of VTK). In particular, VTU allows for
+ * the compression of data and consequently leads to much smaller file
+ * sizes that equivalent VTK files for large files. Since all visualization
+ * programs that support VTK also support VTU, you should consider using the
+ * latter file format instead, by using the write_vtu() function.
*/
template
static void write_vtk (const std::vector > &patches,
@@ -1703,7 +1808,7 @@ public:
/**
- * Write the given list of patches to the output stream in VTK
+ * Write the given list of patches to the output stream in VTU
* format. The data is written in the XML-based VTK format as opposed to the
* traditional format that write_vtk() produces.
*
@@ -1760,7 +1865,35 @@ public:
const std::vector > &vector_data_ranges,
const VtkFlags &flags,
std::ostream &out);
-
+
+ /**
+ * Write the given list of patches to the output stream in SVG format.
+ *
+ * SVG (Scalable Vector Graphics) is an XML-based vector image format
+ * developed and maintained by the World Wide Web Consortium (W3C).
+ * This function conforms to the latest specification SVG 1.1,
+ * released on August 16, 2011. Controlling the graphic output is
+ * possible by setting or clearing the respective flags (see the
+ * SvgFlags struct). At present, this format only supports output
+ * for two-dimensional data, with values in the third direction
+ * taken from a data vector.
+ *
+ * For the output, each patch is subdivided into four triangles
+ * which are then written as polygons and filled with a linear
+ * color gradient. The arising coloring of the patches visualizes
+ * the data values at the vertices taken from the specified data
+ * vector. A colorbar can be drawn to encode the coloring.
+ *
+ * @note Yet only implemented for two dimensions with an additional
+ * dimension reserved for data information.
+ */
+ template
+ static void write_svg (const std::vector > &patches,
+ const std::vector &data_names,
+ const std::vector > &vector_data_ranges,
+ const SvgFlags &flags,
+ std::ostream &out);
+
/**
* Write the given list of patches to the output stream in deal.II
* intermediate format. This is not a format understood by any other
@@ -1909,6 +2042,7 @@ public:
*
tecplot_binary: .plt
*
vtk: .vtk
*
vtu: .vtu
+ *
svg: .svg
*
deal_II_intermediate: .d2.
*
*/
@@ -1988,6 +2122,72 @@ private:
const bool double_precision,
STREAM &out);
+
+ /**
+ * This function projects a three-dimensional point (Point<3> point)
+ * onto a two-dimensional image plane, specified by the position of
+ * the camera viewing system (Point<3> camera_position), camera
+ * direction (Point<3> camera_position), camera horizontal (Point<3>
+ * camera_horizontal, necessary for the correct alignment of the
+ * later images), and the focus of the camera (float camera_focus).
+ *
+ * For SVG output.
+ */
+ static Point<2> svg_project_point(Point<3> point,
+ Point<3> camera_position,
+ Point<3> camera_direction,
+ Point<3> camera_horizontal,
+ float camera_focus);
+ /**
+ * Function to compute the gradient parameters for
+ * a triangle with given values for the vertices.
+ *
+ * Used for svg output.
+ */
+ static Point<6> svg_get_gradient_parameters(Point<3> points[]);
+
+ /**
+ * Class holding the data of one cell of a patch in two space
+ * dimensions for output. It is the projection of a cell in
+ * three-dimensional space (two coordinates, one height value)
+ * to the direction of sight.
+ */
+ class SvgCell
+ {
+ public:
+
+ // Center of the cell (three-dimensional)
+ Point<3> center;
+
+ /**
+ * Vector of vertices of this cell (three-dimensional)
+ */
+ Point<3> vertices[4];
+
+ /**
+ * Depth into the picture, which
+ * is defined as the distance from
+ * an observer at an the origin in
+ * direction of the line of sight.
+ */
+ float depth;
+
+ /**
+ * Vector of vertices of this cell (projected, two-dimensional).
+ */
+ Point<2> projected_vertices[4];
+
+ // Center of the cell (projected, two-dimensional)
+ Point<2> projected_center;
+
+ /**
+ * Comparison operator for
+ * sorting.
+ */
+ bool operator < (const SvgCell &) const;
+ };
+
+
/**
* Class holding the data of one
* cell of a patch in two space
@@ -2290,6 +2490,13 @@ public:
* and write it to out
* in Vtk format. See
* DataOutBase::write_vtk.
+ *
+ * @note VTK is a legacy format and has largely been supplanted by the VTU
+ * format (an XML-structured version of VTK). In particular, VTU allows for
+ * the compression of data and consequently leads to much smaller file
+ * sizes that equivalent VTK files for large files. Since all visualization
+ * programs that support VTK also support VTU, you should consider using the
+ * latter file format instead, by using the write_vtu() function.
*/
void write_vtk (std::ostream &out) const;
@@ -2463,6 +2670,14 @@ public:
void write_visit_record (std::ostream &out,
const std::vector &piece_names) const;
+ /**
+ * Obtain data through get_patches()
+ * and write it to out
+ * in SVG format. See
+ * DataOutBase::write_svg.
+ */
+ void write_svg(std::ostream &out) const;
+
/**
* Obtain data through get_patches()
* and write it to out
@@ -2568,6 +2783,12 @@ public:
*/
void set_flags (const VtkFlags &vtk_flags);
+ /**
+ * Set the flags to be used for
+ * output in SVG format.
+ */
+ void set_flags (const SvgFlags &svg_flags);
+
/**
* Set the flags to be used for output in
* deal.II intermediate format.
@@ -2806,6 +3027,15 @@ private:
*/
VtkFlags vtk_flags;
+ /**
+ * Flags to be used upon output
+ * of svg data in one space
+ * dimension. Can be changed by
+ * using the set_flags
+ * function.
+ */
+ SvgFlags svg_flags;
+
/**
* Flags to be used upon output of
* deal.II intermediate data in one space
diff --git a/deal.II/include/deal.II/base/index_set.h b/deal.II/include/deal.II/base/index_set.h
index c0cf5a1ddd..aff873cba7 100644
--- a/deal.II/include/deal.II/base/index_set.h
+++ b/deal.II/include/deal.II/base/index_set.h
@@ -490,8 +490,8 @@ private:
const Range &range_2)
{
return ((range_1.begin == range_2.begin)
- ||
- (range_1.begin == range_2.begin));
+ &&
+ (range_1.end == range_2.end));
}
std::size_t memory_consumption () const
@@ -571,6 +571,32 @@ private:
};
+/**
+ * Create and return an index set of size $N$ that contains every
+ * single index within this range. In essence, this function
+ * returns an index set created by
+ * @code
+ * IndexSet is (N);
+ * is.add_range(0, N);
+ * @endcode
+ * This function exists so that one can create and initialize
+ * index sets that are complete in one step, or so one can write
+ * code like
+ * @code
+ * if (my_index_set == complete_index_set(my_index_set.size())
+ * ...
+ * @endcode
+ *
+ * @relates IndexSet
+ */
+inline
+IndexSet complete_index_set (const unsigned int N)
+{
+ IndexSet is (N);
+ is.add_range(0, N);
+ return is;
+}
+
/* ------------------ inline functions ------------------ */
inline
diff --git a/deal.II/include/deal.II/base/mpi.h b/deal.II/include/deal.II/base/mpi.h
index 821ad1e4d3..4887776510 100644
--- a/deal.II/include/deal.II/base/mpi.h
+++ b/deal.II/include/deal.II/base/mpi.h
@@ -151,6 +151,14 @@ namespace Utilities
* MPI_Allreduce function, i.e. all processors receive
* the result of this operation.
*
+ * @note Sometimes, not all processors need a results and in that case
+ * one would call the MPI_Reduce function instead of the
+ * MPI_Allreduce function. The latter is at most twice as
+ * expensive, so if you are concerned about performance, it may be
+ * worthwhile investigating whether your algorithm indeed needs the
+ * result everywhere or whether you could get away with calling the
+ * current function and getting the result everywhere.
+ *
* @note This function is only implemented for certain template
* arguments T, namely float, double, int,
* unsigned int.
@@ -198,6 +206,14 @@ namespace Utilities
* MPI_Allreduce function, i.e. all processors receive
* the result of this operation.
*
+ * @note Sometimes, not all processors need a results and in that case
+ * one would call the MPI_Reduce function instead of the
+ * MPI_Allreduce function. The latter is at most twice as
+ * expensive, so if you are concerned about performance, it may be
+ * worthwhile investigating whether your algorithm indeed needs the
+ * result everywhere or whether you could get away with calling the
+ * current function and getting the result everywhere.
+ *
* @note This function is only implemented for certain template
* arguments T, namely float, double, int,
* unsigned int.
@@ -255,12 +271,20 @@ namespace Utilities
* Returns sum, average, minimum,
* maximum, processor id of minimum and
* maximum as a collective operation of
- * on the given MPI communicator @param
+ * on the given MPI communicator @p
* mpi_communicator . Each processor's
- * value is given in @param my_value and
- * the result will be returned in @p
- * result . The result is available on all
+ * value is given in @p my_value and
+ * the result will be returned.
+ * The result is available on all
* machines.
+ *
+ * @note Sometimes, not all processors need a results and in that case
+ * one would call the MPI_Reduce function instead of the
+ * MPI_Allreduce function. The latter is at most twice as
+ * expensive, so if you are concerned about performance, it may be
+ * worthwhile investigating whether your algorithm indeed needs the
+ * result everywhere or whether you could get away with calling the
+ * current function and getting the result everywhere.
*/
MinMaxAvg
min_max_avg (const double my_value,
diff --git a/deal.II/include/deal.II/base/partitioner.h b/deal.II/include/deal.II/base/partitioner.h
index fff826924b..68d4c4d720 100644
--- a/deal.II/include/deal.II/base/partitioner.h
+++ b/deal.II/include/deal.II/base/partitioner.h
@@ -272,7 +272,7 @@ namespace Utilities
* Returns the MPI communicator underlying the
* partitioner object.
*/
- MPI_Comm get_communicator() const;
+ const MPI_Comm& get_communicator() const;
/**
* Computes the memory consumption of this
@@ -565,7 +565,7 @@ namespace Utilities
inline
- MPI_Comm
+ const MPI_Comm&
Partitioner::get_communicator() const
{
return communicator;
diff --git a/deal.II/include/deal.II/base/point.h b/deal.II/include/deal.II/base/point.h
index b91947b89f..c33ea938c8 100644
--- a/deal.II/include/deal.II/base/point.h
+++ b/deal.II/include/deal.II/base/point.h
@@ -53,21 +53,22 @@ class Point : public Tensor<1,dim,Number>
public:
/**
* Standard constructor. Creates
- * an origin.
+ * an object that corresponds to the origin, i.e., all coordinates
+ * are set to zero.
*/
Point ();
/**
* Constructor. Initialize all
* entries to zero if
- * initialize==true.
+ * initialize==true (in which case it is equivalent to the default
+ * constructor) or leaves the elements uninitialized if
+ * initialize==false.
*/
explicit Point (const bool initialize);
/**
- * Convert a tensor to a point. Since no
- * additional data is inside a point,
- * this is ok.
+ * Convert a tensor to a point.
*/
Point (const Tensor<1,dim,Number> &);
@@ -76,7 +77,8 @@ public:
* points. This function is only
* implemented for dim==1 since
* the usage is considered unsafe for
- * points with dim!=1.
+ * points with dim!=1 as it would leave some components
+ * of the point coordinates uninitialized.
*/
explicit Point (const Number x);
@@ -85,18 +87,25 @@ public:
* points. This function is only
* implemented for dim==2 since
* the usage is considered unsafe for
- * points with dim!=2.
+ * points with dim!=2 as it would leave some components
+ * of the point coordinates uninitialized (if dim>2) or would
+ * not use some arguments (if dim<2).
*/
- Point (const Number x, const Number y);
+ Point (const Number x,
+ const Number y);
/**
* Constructor for three dimensional
* points. This function is only
* implemented for dim==3 since
* the usage is considered unsafe for
- * points with dim!=3.
+ * points with dim!=3 as it would leave some components
+ * of the point coordinates uninitialized (if dim>3) or would
+ * not use some arguments (if dim<3).
*/
- Point (const Number x, const Number y, const Number z);
+ Point (const Number x,
+ const Number y,
+ const Number z);
/**
* Return a unit vector in
diff --git a/deal.II/include/deal.II/base/polynomial.h b/deal.II/include/deal.II/base/polynomial.h
index 8897467e04..ec33ec38a8 100644
--- a/deal.II/include/deal.II/base/polynomial.h
+++ b/deal.II/include/deal.II/base/polynomial.h
@@ -30,8 +30,8 @@ DEAL_II_NAMESPACE_OPEN
*/
/**
- * A namespace in which classes relating to the description of
- * 1d polynomial spaces are declared.
+ * A namespace in which classes relating to the description of 1d polynomial
+ * spaces are declared.
*/
namespace Polynomials
{
@@ -50,145 +50,97 @@ namespace Polynomials
{
public:
/**
- * Constructor. The
- * coefficients of the
- * polynomial are passed as
- * arguments, and denote the
- * polynomial $\sum_i a[i]
- * x^i$, i.e. the first element
- * of the array denotes the
- * constant term, the second
- * the linear one, and so
- * on. The degree of the
- * polynomial represented by
- * this object is thus the
- * number of elements in the
- * coefficient array
- * minus one.
+ * Constructor. The coefficients of the polynomial are passed as
+ * arguments, and denote the polynomial $\sum_i a[i] x^i$, i.e. the first
+ * element of the array denotes the constant term, the second the linear
+ * one, and so on. The degree of the polynomial represented by this object
+ * is thus the number of elements in the coefficient array minus
+ * one.
*/
Polynomial (const std::vector &coefficients);
/**
- * Constructor creating a zero
- * polynomial of degree @p n.
+ * Constructor creating a zero polynomial of degree @p n.
*/
Polynomial (const unsigned int n);
/**
- * Constructor for Lagrange polynomial and its
- * point of evaluation. The idea is to
- * construct $\prod_{i\neq j}
- * \frac{x-x_i}{x_j-x_i}$, where j is the
- * evaluation point specified as argument and
- * the support points contain all points
- * (including x_j, which will internally not
- * be stored).
+ * Constructor for Lagrange polynomial and its point of evaluation. The
+ * idea is to construct $\prod_{i\neq j} \frac{x-x_i}{x_j-x_i}$, where j
+ * is the evaluation point specified as argument and the support points
+ * contain all points (including x_j, which will internally not be
+ * stored).
*/
Polynomial (const std::vector > &lagrange_support_points,
const unsigned int evaluation_point);
/**
- * Default constructor creating
- * an illegal object.
+ * Default constructor creating an illegal object.
*/
Polynomial ();
/**
- * Return the value of this
- * polynomial at the given point.
+ * Return the value of this polynomial at the given point.
*
- * This function uses the Horner
- * scheme for numerical stability
- * of the evaluation.
+ * This function uses the Horner scheme for numerical stability of the
+ * evaluation.
*/
number value (const number x) const;
/**
- * Return the values and the
- * derivatives of the
- * Polynomial at point x.
- * values[i],
- * i=0,...,values.size()-1
- * includes the ith
- * derivative. The number of
- * derivatives to be computed is
- * thus determined by the size of
- * the array passed.
+ * Return the values and the derivatives of the Polynomial at point
+ * x. values[i], i=0,...,values.size()-1 includes the
+ * ith derivative. The number of derivatives to be computed is
+ * thus determined by the size of the array passed.
*
- * This function uses the Horner
- * scheme for numerical stability
- * of the evaluation.
+ * This function uses the Horner scheme for numerical stability of the
+ * evaluation.
*/
void value (const number x,
std::vector &values) const;
/**
- * Degree of the polynomial. This
- * is the degree reflected by the
- * number of coefficients
- * provided by the
- * constructor. Leading non-zero
- * coefficients are not treated
- * separately.
+ * Degree of the polynomial. This is the degree reflected by the number of
+ * coefficients provided by the constructor. Leading non-zero coefficients
+ * are not treated separately.
*/
unsigned int degree () const;
/**
- * Scale the abscissa of the
- * polynomial. Given the
- * polynomial p(t) and the
- * scaling t = ax, then the
- * result of this operation is
- * the polynomial q, such that
- * q(x) = p(t).
+ * Scale the abscissa of the polynomial. Given the polynomial p(t)
+ * and the scaling t = ax, then the result of this operation is the
+ * polynomial q, such that q(x) = p(t).
*
- * The operation is performed in
- * place.
+ * The operation is performed in place.
*/
void scale (const number factor);
/**
- * Shift the abscissa oft the
- * polynomial. Given the
- * polynomial p(t) and the
- * shift t = x + a, then the
- * result of this operation is
- * the polynomial q, such that
- * q(x) = p(t).
+ * Shift the abscissa oft the polynomial. Given the polynomial
+ * p(t) and the shift t = x + a, then the result of this
+ * operation is the polynomial q, such that q(x) = p(t).
*
- * The template parameter allows
- * to compute the new
- * coefficients with higher
- * accuracy, since all
- * computations are performed
- * with type number2. This
- * may be necessary, since this
- * operation involves a big
- * number of additions. On a Sun
- * Sparc Ultra with Solaris 2.8,
- * the difference between
- * double and long double
- * was not significant, though.
+ * The template parameter allows to compute the new coefficients with
+ * higher accuracy, since all computations are performed with type
+ * number2. This may be necessary, since this operation involves
+ * a big number of additions. On a Sun Sparc Ultra with Solaris 2.8, the
+ * difference between double and long double was not
+ * significant, though.
*
- * The operation is performed in
- * place, i.e. the coefficients
- * of the present object are
- * changed.
+ * The operation is performed in place, i.e. the coefficients of the
+ * present object are changed.
*/
template
void shift (const number2 offset);
/**
- * Compute the derivative of a
- * polynomial.
+ * Compute the derivative of a polynomial.
*/
Polynomial derivative () const;
/**
- * Compute the primitive of a
- * polynomial. the coefficient
- * of the zero order term of
- * the polynomial is zero.
+ * Compute the primitive of a polynomial. the coefficient of the zero
+ * order term of the polynomial is zero.
*/
Polynomial primitive () const;
@@ -213,8 +165,8 @@ namespace Polynomials
Polynomial &operator -= (const Polynomial &p);
/**
- * Test for equality of two polynomials.
- */
+ * Test for equality of two polynomials.
+ */
bool operator == (const Polynomial &p) const;
/**
@@ -223,8 +175,8 @@ namespace Polynomials
void print(std::ostream &out) const;
/**
- * Write or read the data of this object to or
- * from a stream for the purpose of serialization.
+ * Write or read the data of this object to or from a stream for the
+ * purpose of serialization.
*/
template
void serialize (Archive &ar, const unsigned int version);
@@ -232,15 +184,13 @@ namespace Polynomials
protected:
/**
- * This function performs the
- * actual scaling.
+ * This function performs the actual scaling.
*/
static void scale (std::vector &coefficients,
const number factor);
/**
- * This function performs the
- * actual shift
+ * This function performs the actual shift
*/
template
static void shift (std::vector &coefficients,
@@ -253,49 +203,37 @@ namespace Polynomials
const number factor);
/**
- * Transforms polynomial form of
- * product of linear factors into
- * standard form, $\sum_i a_i
- * x^i$. Deletes all data structures
- * related to the product form.
+ * Transforms polynomial form of product of linear factors into standard
+ * form, $\sum_i a_i x^i$. Deletes all data structures related to the
+ * product form.
*/
void transform_into_standard_form ();
/**
- * Coefficients of the polynomial
- * $\sum_i a_i x^i$. This vector
- * is filled by the constructor
- * of this class and may be
- * passed down by derived
+ * Coefficients of the polynomial $\sum_i a_i x^i$. This vector is filled
+ * by the constructor of this class and may be passed down by derived
* classes.
*
- * This vector cannot be constant
- * since we want to allow copying
- * of polynomials.
+ * This vector cannot be constant since we want to allow copying of
+ * polynomials.
*/
std::vector coefficients;
/**
- * Stores whether the polynomial is in
- * Lagrange product form, i.e.,
- * constructed as a product $(x-x_0)
- * (x-x_1) \ldots (x-x_n)/c$, or not.
+ * Stores whether the polynomial is in Lagrange product form, i.e.,
+ * constructed as a product $(x-x_0) (x-x_1) \ldots (x-x_n)/c$, or not.
*/
bool in_lagrange_product_form;
/**
- * If the polynomial is in Lagrange
- * product form, i.e., constructed as a
- * product $(x-x_0) (x-x_1) \ldots
- * (x-x_n)/c$, store the shifts $x_i$.
+ * If the polynomial is in Lagrange product form, i.e., constructed as a
+ * product $(x-x_0) (x-x_1) \ldots (x-x_n)/c$, store the shifts $x_i$.
*/
std::vector lagrange_support_points;
/**
- * If the polynomial is in Lagrange
- * product form, i.e., constructed as a
- * product $(x-x_0) (x-x_1) \ldots
- * (x-x_n)/c$, store the weight c.
+ * If the polynomial is in Lagrange product form, i.e., constructed as a
+ * product $(x-x_0) (x-x_1) \ldots (x-x_n)/c$, store the weight c.
*/
number lagrange_weight;
};
@@ -367,32 +305,20 @@ namespace Polynomials
{
public:
/**
- * Constructor. Takes the degree
- * n of the Lagrangian
- * polynom and the index
- * support_point of the
- * support point. Fills the
- * coefficients of the base
- * class Polynomial.
+ * Constructor. Takes the degree n of the Lagrangian polynom and
+ * the index support_point of the support point. Fills the
+ * coefficients of the base class Polynomial.
*/
LagrangeEquidistant (const unsigned int n,
const unsigned int support_point);
/**
- * Return a vector of polynomial
- * objects of degree degree,
- * which then spans the full
- * space of polynomials up to the
- * given degree. The polynomials
- * are generated by calling the
- * constructor of this class with
- * the same degree but support
- * point running from zero to
- * degree. This function may
- * be used to initialize the
- * TensorProductPolynomials
- * and PolynomialSpace
- * classes.
+ * Return a vector of polynomial objects of degree degree, which
+ * then spans the full space of polynomials up to the given degree. The
+ * polynomials are generated by calling the constructor of this class with
+ * the same degree but support point running from zero to
+ * degree. This function may be used to initialize the
+ * TensorProductPolynomials and PolynomialSpace classes.
*/
static
std::vector >
@@ -401,12 +327,8 @@ namespace Polynomials
private:
/**
- * Computes the coefficients
- * of the base class
- * Polynomial. This function
- * is static to allow to be
- * called in the
- * constructor.
+ * Computes the coefficients of the base class Polynomial. This
+ * function is static to allow to be called in the constructor.
*/
static
void
@@ -416,15 +338,12 @@ namespace Polynomials
};
+
/**
- * Given a set of points along the
- * real axis, this function returns
- * all Lagrange polynomials for
- * interpolation of these
- * points. The number of
- * polynomials is equal to the
- * number of points and the maximum
- * degree is one less.
+ * Given a set of points along the real axis, this function returns all
+ * Lagrange polynomials for interpolation of these points. The number of
+ * polynomials is equal to the number of points and the maximum degree is
+ * one less.
*/
std::vector >
generate_complete_Lagrange_basis (const std::vector > &points);
@@ -801,8 +720,7 @@ namespace Polynomials
void
Polynomial::serialize (Archive &ar, const unsigned int)
{
- // forward to serialization
- // function in the base class.
+ // forward to serialization function in the base class.
ar &static_cast(*this);
ar &coefficients;
ar &in_lagrange_product_form;
diff --git a/deal.II/include/deal.II/base/polynomials_piecewise.h b/deal.II/include/deal.II/base/polynomials_piecewise.h
new file mode 100644
index 0000000000..c9e1dbbcaa
--- /dev/null
+++ b/deal.II/include/deal.II/base/polynomials_piecewise.h
@@ -0,0 +1,224 @@
+//---------------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009, 2010, 2011, 2012 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__polynomials_piecewise_h
+#define __deal2__polynomials_piecewise_h
+
+
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * @addtogroup Polynomials
+ * @{
+ */
+
+/**
+ * A namespace in which classes relating to the description of 1d polynomial
+ * spaces are declared.
+ */
+namespace Polynomials
+{
+
+ /**
+ * Definition of piecewise 1D polynomials for the unit interval. This space
+ * allows the description of interpolating polynomials on parts of the unit
+ * interval, similarly to the definition of finite element basis functions
+ * on the subdivided elements. This primary purpose of this class is to
+ * allow constructing FE_Q_iso_Q1 elements that put additional degrees of
+ * freedom into an equivalent of a refined mesh instead of higher order
+ * polynomials, which is useful when using mixed finite elements.
+ *
+ * @author Martin Kronbichler, 2013
+ */
+ template
+ class PiecewisePolynomial : public Subscriptor
+ {
+ public:
+ /**
+ * Constructor for Lagrange polynomial on an interval that is a subset of
+ * the unit interval. It uses a polynomial description that is scaled to
+ * the size of the subinterval compared to the unit interval, the total
+ * number of intervals (subdivisions), the current index of the interval
+ * as well as if the polynomial spans onto the next interval (e.g., if it
+ * lives on two neighboring intervals).
+ *
+ * If the number of intervals is one, the piecewise polynomial behaves
+ * exactly like a usual polynomial.
+ */
+ PiecewisePolynomial (const Polynomial &coefficients_on_interval,
+ const unsigned int n_intervals,
+ const unsigned int interval,
+ const bool spans_next_interval);
+
+ /**
+ * Return the value of this polynomial at the given point, evaluating the
+ * underlying polynomial. The polynomial evaluates to zero when outside of
+ * the given interval (and possible the next one to the right when it
+ * spans over that range).
+ */
+ number value (const number x) const;
+
+ /**
+ * Return the values and the derivatives of the Polynomial at point
+ * x. values[i], i=0,...,values.size()-1 includes the
+ * ith derivative. The number of derivatives to be computed is
+ * thus determined by the size of the array passed.
+ *
+ * Note that all the derivatives evaluate to zero at the border between
+ * intervals (assuming exact arithmetics) in the interior of the unit
+ * interval, as there is no unique gradient value in that case for a
+ * piecewise polynomial. This is not always desired (e.g., when evaluating
+ * jumps of gradients on the element boundary), but it is the user's
+ * responsibility to avoid evaluation at these points when it does not
+ * make sense.
+ */
+ void value (const number x,
+ std::vector &values) const;
+
+ /**
+ * Degree of the polynomial. This is the degree of the underlying base
+ * polynomial.
+ */
+ unsigned int degree () const;
+
+ /**
+ * Write or read the data of this object to or from a stream for the
+ * purpose of serialization.
+ */
+ template
+ void serialize (Archive &ar, const unsigned int version);
+
+ protected:
+
+ /**
+ * Underlying polynomial object that is scaled to a subinterval and
+ * concatenated accordingly.
+ */
+ Polynomial polynomial;
+
+ /**
+ * Stores the number of intervals that the unit interval is divided into.
+ */
+ unsigned int n_intervals;
+
+ /**
+ * Stores the index of the current polynomial in the range of
+ * intervals.
+ */
+ unsigned int interval;
+
+ /**
+ * Store if the polynomial spans over two adjacent intervals, i.e., the
+ * one given in subinterval and the next one.
+ */
+ bool spans_two_intervals;
+ };
+
+
+
+ /**
+ * Generates a complete Lagrange basis on a subdivision of the unit interval
+ * in smaller intervals for a given degree on the subintervals and number of
+ * intervals.
+ */
+ std::vector >
+ generate_complete_Lagrange_basis_on_subdivisions (const unsigned int n_subdivisions,
+ const unsigned int base_degree);
+
+}
+
+
+/** @} */
+
+/* -------------------------- inline functions --------------------- */
+
+namespace Polynomials
+{
+ template
+ inline
+ unsigned int
+ PiecewisePolynomial::degree () const
+ {
+ return polynomial.degree();
+ }
+
+
+
+ template
+ inline
+ number
+ PiecewisePolynomial::value (const number x) const
+ {
+ AssertIndexRange (interval, n_intervals);
+ number y = x;
+ // shift polynomial if necessary
+ if (n_intervals > 1)
+ {
+ const number step = 1./n_intervals;
+
+ // polynomial spans over two intervals
+ if (spans_two_intervals == true)
+ {
+ const number offset = step * interval;
+ if (xoffset+step+step)
+ return 0;
+ else if (xoffset+step)
+ return 0;
+ else
+ y = x-offset;
+ }
+
+ return polynomial.value(y);
+ }
+ else
+ return polynomial.value(x);
+ }
+
+
+
+ template
+ template
+ inline
+ void
+ PiecewisePolynomial::serialize (Archive &ar, const unsigned int)
+ {
+ // forward to serialization function in the base class.
+ ar &static_cast(*this);
+ ar &polynomial;
+ ar &n_intervals;
+ ar &interval;
+ ar &spans_two_intervals;
+ }
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/include/deal.II/base/tensor_product_polynomials.h b/deal.II/include/deal.II/base/tensor_product_polynomials.h
index 972f0ab0dd..3babe00729 100644
--- a/deal.II/include/deal.II/base/tensor_product_polynomials.h
+++ b/deal.II/include/deal.II/base/tensor_product_polynomials.h
@@ -58,80 +58,57 @@ DEAL_II_NAMESPACE_OPEN
*
* @author Ralf Hartmann, 2000, 2004, Guido Kanschat, 2000, Wolfgang Bangerth 2003
*/
-template
+template >
class TensorProductPolynomials
{
public:
/**
- * Access to the dimension of
- * this object, for checking and
- * automatic setting of dimension
- * in other classes.
+ * Access to the dimension of this object, for checking and automatic
+ * setting of dimension in other classes.
*/
static const unsigned int dimension = dim;
/**
- * Constructor. pols is
- * a vector of objects that
- * should be derived or otherwise
- * convertible to one-dimensional
- * polynomial objects. It will be
- * copied element by element into
- * a private variable.
+ * Constructor. pols is a vector of objects that should be derived
+ * or otherwise convertible to one-dimensional polynomial objects of type @p
+ * POLY (template argument of class). It will be copied element by element
+ * into a private variable.
*/
template
TensorProductPolynomials (const std::vector &pols);
/**
- * Prints the list of the indices
- * to out.
+ * Prints the list of the indices to out.
*/
void output_indices(std::ostream &out) const;
/**
- * Sets the ordering of the
- * polynomials. Requires
- * renumber.size()==n().
- * Stores a copy of
- * renumber.
+ * Sets the ordering of the polynomials. Requires
+ * renumber.size()==n(). Stores a copy of renumber.
*/
void set_numbering(const std::vector &renumber);
/**
- * Gives read access to the
- * renumber vector.
+ * Gives read access to the renumber vector.
*/
const std::vector &get_numbering() const;
/**
- * Gives read access to the
- * inverse renumber vector.
+ * Gives read access to the inverse renumber vector.
*/
const std::vector &get_numbering_inverse() const;
/**
- * Computes the value and the
- * first and second derivatives
- * of each tensor product
- * polynomial at unit_point.
+ * Computes the value and the first and second derivatives of each tensor
+ * product polynomial at unit_point.
*
- * The size of the vectors must
- * either be equal 0 or equal
- * n(). In the first case, the
- * function will not compute
- * these values.
+ * The size of the vectors must either be equal 0 or equal n(). In the first
+ * case, the function will not compute these values.
*
- * If you need values or
- * derivatives of all tensor
- * product polynomials then use
- * this function, rather than
- * using any of the
- * compute_value(),
- * compute_grad() or
- * compute_grad_grad()
- * functions, see below, in a
- * loop over all tensor product
- * polynomials.
+ * If you need values or derivatives of all tensor product polynomials then
+ * use this function, rather than using any of the compute_value(),
+ * compute_grad() or compute_grad_grad() functions, see below, in a loop
+ * over all tensor product polynomials.
*/
void compute (const Point &unit_point,
std::vector &values,
@@ -139,157 +116,89 @@ public:
std::vector > &grad_grads) const;
/**
- * Computes the value of the
- * ith tensor product
- * polynomial at
- * unit_point. Here i is
- * given in tensor product
+ * Computes the value of the ith tensor product polynomial at
+ * unit_point. Here i is given in tensor product
* numbering.
*
- * Note, that using this function
- * within a loop over all tensor
- * product polynomials is not
- * efficient, because then each
- * point value of the underlying
- * (one-dimensional) polynomials
- * is (unnecessarily) computed
- * several times. Instead use
- * the compute() function with
- * values.size()==n()
- * to get the point values of all
- * tensor polynomials all at once
- * and in a much more efficient
- * way.
+ * Note, that using this function within a loop over all tensor product
+ * polynomials is not efficient, because then each point value of the
+ * underlying (one-dimensional) polynomials is (unnecessarily) computed
+ * several times. Instead use the compute() function with
+ * values.size()==n() to get the point values of all tensor
+ * polynomials all at once and in a much more efficient way.
*/
double compute_value (const unsigned int i,
const Point &p) const;
/**
- * Computes the grad of the
- * ith tensor product
- * polynomial at
- * unit_point. Here i is
- * given in tensor product
+ * Computes the grad of the ith tensor product polynomial at
+ * unit_point. Here i is given in tensor product
* numbering.
*
- * Note, that using this function
- * within a loop over all tensor
- * product polynomials is not
- * efficient, because then each
- * derivative value of the
- * underlying (one-dimensional)
- * polynomials is (unnecessarily)
- * computed several times.
- * Instead use the compute()
- * function, see above, with
- * grads.size()==n()
- * to get the point value of all
- * tensor polynomials all at once
- * and in a much more efficient
- * way.
+ * Note, that using this function within a loop over all tensor product
+ * polynomials is not efficient, because then each derivative value of the
+ * underlying (one-dimensional) polynomials is (unnecessarily) computed
+ * several times. Instead use the compute() function, see above, with
+ * grads.size()==n() to get the point value of all tensor
+ * polynomials all at once and in a much more efficient way.
*/
Tensor<1,dim> compute_grad (const unsigned int i,
const Point &p) const;
/**
- * Computes the second
- * derivative (grad_grad) of the
- * ith tensor product
- * polynomial at
- * unit_point. Here i is
- * given in tensor product
- * numbering.
+ * Computes the second derivative (grad_grad) of the ith tensor
+ * product polynomial at unit_point. Here i is given in
+ * tensor product numbering.
*
- * Note, that using this function
- * within a loop over all tensor
- * product polynomials is not
- * efficient, because then each
- * derivative value of the
- * underlying (one-dimensional)
- * polynomials is (unnecessarily)
- * computed several times.
- * Instead use the compute()
- * function, see above, with
- * grad_grads.size()==n()
- * to get the point value of all
- * tensor polynomials all at once
- * and in a much more efficient
- * way.
+ * Note, that using this function within a loop over all tensor product
+ * polynomials is not efficient, because then each derivative value of the
+ * underlying (one-dimensional) polynomials is (unnecessarily) computed
+ * several times. Instead use the compute() function, see above, with
+ * grad_grads.size()==n() to get the point value of all tensor
+ * polynomials all at once and in a much more efficient way.
*/
Tensor<2,dim> compute_grad_grad (const unsigned int i,
const Point &p) const;
/**
- * Returns the number of tensor
- * product polynomials. For n
- * 1d polynomials this is ndim.
+ * Returns the number of tensor product polynomials. For n 1d
+ * polynomials this is ndim.
*/
unsigned int n () const;
-private:
+protected:
/**
- * Copy of the vector pols of
- * polynomials given to the
- * constructor.
+ * Copy of the vector pols of polynomials given to the constructor.
*/
- std::vector > polynomials;
+ std::vector polynomials;
/**
- * Number of tensor product
- * polynomials. See n().
+ * Number of tensor product polynomials. See n().
*/
unsigned int n_tensor_pols;
/**
- * Index map for reordering the
- * polynomials.
+ * Index map for reordering the polynomials.
*/
std::vector index_map;
/**
- * Index map for reordering the
- * polynomials.
+ * Index map for reordering the polynomials.
*/
std::vector index_map_inverse;
/**
- * Each tensor product polynomial
- * i is a product of
- * one-dimensional polynomials in
- * each space direction. Compute
- * the indices of these
- * one-dimensional polynomials
- * for each space direction,
- * given the index i.
+ * Each tensor product polynomial i is a product of one-dimensional
+ * polynomials in each space direction. Compute the indices of these
+ * one-dimensional polynomials for each space direction, given the index
+ * i.
*/
- // fix to avoid compiler warnings about zero
- // length arrays
+ // fix to avoid compiler warnings about zero length arrays
void compute_index (const unsigned int i,
unsigned int (&indices)[(dim>0?dim:1)]) const;
};
-#ifndef DOXYGEN
-
-template
-inline
-const std::vector &
-TensorProductPolynomials::get_numbering() const
-{
- return index_map;
-}
-
-
-template
-inline
-const std::vector &
-TensorProductPolynomials::get_numbering_inverse() const
-{
- return index_map_inverse;
-}
-
-
-#endif // DOXYGEN
/**
@@ -316,44 +225,26 @@ class AnisotropicPolynomials
{
public:
/**
- * Constructor. pols is a
- * table of one-dimensional
- * polynomials. The number of
- * rows in this table should be
- * equal to the space dimension,
- * with the elements of each row
- * giving the polynomials that
- * shall be used in this
- * particular coordinate
- * direction. These polynomials
- * may vary between coordinates,
- * as well as their number.
+ * Constructor. pols is a table of one-dimensional polynomials. The
+ * number of rows in this table should be equal to the space dimension, with
+ * the elements of each row giving the polynomials that shall be used in
+ * this particular coordinate direction. These polynomials may vary between
+ * coordinates, as well as their number.
*/
AnisotropicPolynomials (const std::vector > > &pols);
/**
- * Computes the value and the
- * first and second derivatives
- * of each tensor product
- * polynomial at unit_point.
+ * Computes the value and the first and second derivatives of each tensor
+ * product polynomial at unit_point.
*
- * The size of the vectors must
- * either be equal 0 or equal
- * n_tensor_pols. In the
- * first case, the function will
- * not compute these values.
+ * The size of the vectors must either be equal 0 or equal
+ * n_tensor_pols. In the first case, the function will not compute
+ * these values.
*
- * If you need values or
- * derivatives of all tensor
- * product polynomials then use
- * this function, rather than
- * using any of the
- * compute_value,
- * compute_grad or
- * compute_grad_grad
- * functions, see below, in a
- * loop over all tensor product
- * polynomials.
+ * If you need values or derivatives of all tensor product polynomials then
+ * use this function, rather than using any of the compute_value,
+ * compute_grad or compute_grad_grad functions, see below,
+ * in a loop over all tensor product polynomials.
*/
void compute (const Point &unit_point,
std::vector &values,
@@ -361,130 +252,79 @@ public:
std::vector > &grad_grads) const;
/**
- * Computes the value of the
- * ith tensor product
- * polynomial at
- * unit_point. Here i is
- * given in tensor product
+ * Computes the value of the ith tensor product polynomial at
+ * unit_point. Here i is given in tensor product
* numbering.
*
- * Note, that using this function
- * within a loop over all tensor
- * product polynomials is not
- * efficient, because then each
- * point value of the underlying
- * (one-dimensional) polynomials
- * is (unnecessarily) computed
- * several times. Instead use
- * the compute function, see
- * above, with
- * values.size()==n_tensor_pols
- * to get the point values of all
- * tensor polynomials all at once
- * and in a much more efficient
- * way.
+ * Note, that using this function within a loop over all tensor product
+ * polynomials is not efficient, because then each point value of the
+ * underlying (one-dimensional) polynomials is (unnecessarily) computed
+ * several times. Instead use the compute function, see above,
+ * with values.size()==n_tensor_pols to get the point values of all
+ * tensor polynomials all at once and in a much more efficient way.
*/
double compute_value (const unsigned int i,
const Point &p) const;
/**
- * Computes the grad of the
- * ith tensor product
- * polynomial at
- * unit_point. Here i is
- * given in tensor product
+ * Computes the grad of the ith tensor product polynomial at
+ * unit_point. Here i is given in tensor product
* numbering.
*
- * Note, that using this function
- * within a loop over all tensor
- * product polynomials is not
- * efficient, because then each
- * derivative value of the
- * underlying (one-dimensional)
- * polynomials is (unnecessarily)
- * computed several times.
- * Instead use the compute
- * function, see above, with
- * grads.size()==n_tensor_pols
- * to get the point value of all
- * tensor polynomials all at once
- * and in a much more efficient
- * way.
+ * Note, that using this function within a loop over all tensor product
+ * polynomials is not efficient, because then each derivative value of the
+ * underlying (one-dimensional) polynomials is (unnecessarily) computed
+ * several times. Instead use the compute function, see above,
+ * with grads.size()==n_tensor_pols to get the point value of all
+ * tensor polynomials all at once and in a much more efficient way.
*/
Tensor<1,dim> compute_grad (const unsigned int i,
const Point &p) const;
/**
- * Computes the second
- * derivative (grad_grad) of the
- * ith tensor product
- * polynomial at
- * unit_point. Here i is
- * given in tensor product
- * numbering.
+ * Computes the second derivative (grad_grad) of the ith tensor
+ * product polynomial at unit_point. Here i is given in
+ * tensor product numbering.
*
- * Note, that using this function
- * within a loop over all tensor
- * product polynomials is not
- * efficient, because then each
- * derivative value of the
- * underlying (one-dimensional)
- * polynomials is (unnecessarily)
- * computed several times.
- * Instead use the compute
- * function, see above, with
- * grad_grads.size()==n_tensor_pols
- * to get the point value of all
- * tensor polynomials all at once
- * and in a much more efficient
- * way.
+ * Note, that using this function within a loop over all tensor product
+ * polynomials is not efficient, because then each derivative value of the
+ * underlying (one-dimensional) polynomials is (unnecessarily) computed
+ * several times. Instead use the compute function, see above,
+ * with grad_grads.size()==n_tensor_pols to get the point value of
+ * all tensor polynomials all at once and in a much more efficient way.
*/
Tensor<2,dim> compute_grad_grad (const unsigned int i,
const Point &p) const;
/**
- * Returns the number of tensor
- * product polynomials. It is the
- * product of the number of
- * polynomials in each coordinate
- * direction.
+ * Returns the number of tensor product polynomials. It is the product of
+ * the number of polynomials in each coordinate direction.
*/
unsigned int n () const;
private:
/**
- * Copy of the vector pols of
- * polynomials given to the
- * constructor.
+ * Copy of the vector pols of polynomials given to the constructor.
*/
std::vector > > polynomials;
/**
- * Number of tensor product
- * polynomials. This is
- * Nx*Ny*Nz, or with terms
- * dropped if the number of space
- * dimensions is less than 3.
+ * Number of tensor product polynomials. This is Nx*Ny*Nz, or with
+ * terms dropped if the number of space dimensions is less than 3.
*/
unsigned int n_tensor_pols;
/**
- * Each tensor product polynomial
- * @þ{i} is a product of
- * one-dimensional polynomials in
- * each space direction. Compute
- * the indices of these
- * one-dimensional polynomials
- * for each space direction,
- * given the index i.
+ * Each tensor product polynomial @þ{i} is a product of one-dimensional
+ * polynomials in each space direction. Compute the indices of these
+ * one-dimensional polynomials for each space direction, given the index
+ * i.
*/
void compute_index (const unsigned int i,
unsigned int (&indices)[dim]) const;
/**
- * Given the input to the
- * constructor, compute
- * n_tensor_pols.
+ * Given the input to the constructor, compute n_tensor_pols.
*/
static
unsigned int
@@ -495,28 +335,14 @@ private:
#ifndef DOXYGEN
-/* -------------- declaration of explicit specializations --- */
-
-template <>
-void
-TensorProductPolynomials<1>::compute_index(const unsigned int n,
- unsigned int (&index)[1]) const;
-template <>
-void
-TensorProductPolynomials<2>::compute_index(const unsigned int n,
- unsigned int (&index)[2]) const;
-template <>
-void
-TensorProductPolynomials<3>::compute_index(const unsigned int n,
- unsigned int (&index)[3]) const;
-
/* ---------------- template and inline functions ---------- */
-template
+template
template
-TensorProductPolynomials::
+inline
+TensorProductPolynomials::
TensorProductPolynomials(const std::vector &pols)
:
polynomials (pols.begin(), pols.end()),
@@ -524,10 +350,8 @@ TensorProductPolynomials(const std::vector &pols)
index_map(n_tensor_pols),
index_map_inverse(n_tensor_pols)
{
- // per default set this index map
- // to identity. This map can be
- // changed by the user through the
- // set_numbering() function
+ // per default set this index map to identity. This map can be changed by
+ // the user through the set_numbering() function
for (unsigned int i=0; i &pols)
-template <>
-void
-AnisotropicPolynomials<1>::compute_index(const unsigned int n,
- unsigned int (&index)[1]) const;
-template <>
-void
-AnisotropicPolynomials<2>::compute_index(const unsigned int n,
- unsigned int (&index)[2]) const;
-template <>
-void
-AnisotropicPolynomials<3>::compute_index(const unsigned int n,
- unsigned int (&index)[3]) const;
+template
+inline
+unsigned int
+TensorProductPolynomials::n() const
+{
+ if (dim == 0)
+ return numbers::invalid_unsigned_int;
+ else
+ return n_tensor_pols;
+}
+
+
+
+template
+inline
+const std::vector &
+TensorProductPolynomials::get_numbering() const
+{
+ return index_map;
+}
+
+
+template
+inline
+const std::vector &
+TensorProductPolynomials::get_numbering_inverse() const
+{
+ return index_map_inverse;
+}
+
+
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/include/deal.II/base/tensor_product_polynomials_const.h b/deal.II/include/deal.II/base/tensor_product_polynomials_const.h
index 9d47b28d72..82fa269e5f 100644
--- a/deal.II/include/deal.II/base/tensor_product_polynomials_const.h
+++ b/deal.II/include/deal.II/base/tensor_product_polynomials_const.h
@@ -18,6 +18,7 @@
#include
#include
#include
+#include
#include
#include
@@ -25,45 +26,21 @@
DEAL_II_NAMESPACE_OPEN
-//TODO: This class is an almost exact copy of the TensorProductPolynomials class. They should be unified. A simple way to do that would be to derive this one from the other and in the implemented member functions simply defer to the corresponding function in the base class whenever appropriate.
-
/**
* @addtogroup Polynomials
* @{
*/
/**
- * Tensor product of given polynomials.
- *
- * Given a vector of n one-dimensional polynomials
- * P1 to Pn, this class generates
- * ndim polynomials of the form
- * Qijk(x,y,z) =
- * Pi(x)Pj(y)Pk(z)and a locally
- * constant function. If the base polynomials are mutually orthogonal
- * on the interval [-1,1] or [0,1], then the tensor product
- * polynomials are orthogonal on [-1,1]dim or [0,1]
- * dim, respectively.
- *
- * Indexing is as follows: the order of dim-dimensional polynomials is
- * x-coordinates running fastest, then y-coordinate, etc. The first
- * few polynomials are thus P1(x)P1(y),
- * P2(x)P1(y), P3(x)P1(y),
- * ..., P1(x)P2(y),
- * P2(x)P2(y), P3(x)P2(y),
- * ... and likewise in 3d. The locally constant function has the last index
- *
- * The output_indices() function prints the ordering of the
- * dim-dimensional polynomials, i.e. for each polynomial in the
- * polynomial space it gives the indices i,j,k of the one-dimensional
- * polynomials in x,y and z direction. The ordering of the
- * dim-dimensional polynomials can be changed by using the
- * set_numbering() function.
+ * Tensor product of given polynomials and a locally constant function. This
+ * class inherits most of its functionality from TensorProductPolynomials. It
+ * works similarly to that class but adds a constant function for the last
+ * index.
*
- * @author Ralf Hartmann, 2000, 2004, Guido Kanschat, 2000, Wolfgang Bangerth 2003
+ * @author Timo Heister, 2012
*/
template
-class TensorProductPolynomialsConst
+class TensorProductPolynomialsConst : public TensorProductPolynomials
{
public:
/**
@@ -75,67 +52,24 @@ public:
static const unsigned int dimension = dim;
/**
- * Constructor. pols is
- * a vector of objects that
- * should be derived or otherwise
- * convertible to one-dimensional
- * polynomial objects. It will be
- * copied element by element into
- * a private variable.
+ * Constructor. pols is a vector of objects that should be derived
+ * or otherwise convertible to one-dimensional polynomial objects. It will
+ * be copied element by element into a private variable.
*/
template
TensorProductPolynomialsConst (const std::vector &pols);
/**
- * Prints the list of the indices
- * to out.
- */
- void output_indices(std::ostream &out) const;
-
- /**
- * Sets the ordering of the
- * polynomials. Requires
- * renumber.size()==n().
- * Stores a copy of
- * renumber.
- */
- void set_numbering(const std::vector &renumber);
-
- /**
- * Gives read access to the
- * renumber vector.
- */
- const std::vector &get_numbering() const;
-
- /**
- * Gives read access to the
- * inverse renumber vector.
- */
- const std::vector &get_numbering_inverse() const;
-
- /**
- * Computes the value and the
- * first and second derivatives
- * of each tensor product
- * polynomial at unit_point.
+ * Computes the value and the first and second derivatives of each tensor
+ * product polynomial at unit_point.
*
- * The size of the vectors must
- * either be equal 0 or equal
- * n(). In the first case, the
- * function will not compute
- * these values.
+ * The size of the vectors must either be equal 0 or equal n(). In the first
+ * case, the function will not compute these values.
*
- * If you need values or
- * derivatives of all tensor
- * product polynomials then use
- * this function, rather than
- * using any of the
- * compute_value(),
- * compute_grad() or
- * compute_grad_grad()
- * functions, see below, in a
- * loop over all tensor product
- * polynomials.
+ * If you need values or derivatives of all tensor product polynomials then
+ * use this function, rather than using any of the compute_value(),
+ * compute_grad() or compute_grad_grad() functions, see below, in a loop
+ * over all tensor product polynomials.
*/
void compute (const Point &unit_point,
std::vector &values,
@@ -143,202 +77,98 @@ public:
std::vector > &grad_grads) const;
/**
- * Computes the value of the
- * ith tensor product
- * polynomial at
- * unit_point. Here i is
- * given in tensor product
+ * Computes the value of the ith tensor product polynomial at
+ * unit_point. Here i is given in tensor product
* numbering.
*
- * Note, that using this function
- * within a loop over all tensor
- * product polynomials is not
- * efficient, because then each
- * point value of the underlying
- * (one-dimensional) polynomials
- * is (unnecessarily) computed
- * several times. Instead use
- * the compute() function with
- * values.size()==n()
- * to get the point values of all
- * tensor polynomials all at once
- * and in a much more efficient
- * way.
+ * Note, that using this function within a loop over all tensor product
+ * polynomials is not efficient, because then each point value of the
+ * underlying (one-dimensional) polynomials is (unnecessarily) computed
+ * several times. Instead use the compute() function with
+ * values.size()==n() to get the point values of all tensor
+ * polynomials all at once and in a much more efficient way.
*/
double compute_value (const unsigned int i,
const Point &p) const;
/**
- * Computes the grad of the
- * ith tensor product
- * polynomial at
- * unit_point. Here i is
- * given in tensor product
+ * Computes the grad of the ith tensor product polynomial at
+ * unit_point. Here i is given in tensor product
* numbering.
*
- * Note, that using this function
- * within a loop over all tensor
- * product polynomials is not
- * efficient, because then each
- * derivative value of the
- * underlying (one-dimensional)
- * polynomials is (unnecessarily)
- * computed several times.
- * Instead use the compute()
- * function, see above, with
- * grads.size()==n()
- * to get the point value of all
- * tensor polynomials all at once
- * and in a much more efficient
- * way.
+ * Note, that using this function within a loop over all tensor product
+ * polynomials is not efficient, because then each derivative value of the
+ * underlying (one-dimensional) polynomials is (unnecessarily) computed
+ * several times. Instead use the compute() function, see above, with
+ * grads.size()==n() to get the point value of all tensor
+ * polynomials all at once and in a much more efficient way.
*/
Tensor<1,dim> compute_grad (const unsigned int i,
const Point &p) const;
/**
- * Computes the second
- * derivative (grad_grad) of the
- * ith tensor product
- * polynomial at
- * unit_point. Here i is
- * given in tensor product
- * numbering.
+ * Computes the second derivative (grad_grad) of the ith tensor
+ * product polynomial at unit_point. Here i is given in
+ * tensor product numbering.
*
- * Note, that using this function
- * within a loop over all tensor
- * product polynomials is not
- * efficient, because then each
- * derivative value of the
- * underlying (one-dimensional)
- * polynomials is (unnecessarily)
- * computed several times.
- * Instead use the compute()
- * function, see above, with
- * grad_grads.size()==n()
- * to get the point value of all
- * tensor polynomials all at once
- * and in a much more efficient
- * way.
+ * Note, that using this function within a loop over all tensor product
+ * polynomials is not efficient, because then each derivative value of the
+ * underlying (one-dimensional) polynomials is (unnecessarily) computed
+ * several times. Instead use the compute() function, see above, with
+ * grad_grads.size()==n() to get the point value of all tensor
+ * polynomials all at once and in a much more efficient way.
*/
Tensor<2,dim> compute_grad_grad (const unsigned int i,
const Point &p) const;
/**
- * Returns the number of tensor
- * product polynomials plus the constant
- * function. For n 1d polynomials
- * this is ndim+1.
+ * Returns the number of tensor product polynomials plus the constant
+ * function. For n 1d polynomials this is ndim+1.
*/
unsigned int n () const;
-
-
-private:
- /**
- * Copy of the vector pols of
- * polynomials given to the
- * constructor.
- */
- std::vector > polynomials;
-
- /**
- * Number of tensor product
- * polynomials. See n().
- */
- unsigned int n_tensor_pols;
-
- /**
- * Index map for reordering the
- * polynomials.
- */
- std::vector index_map;
-
- /**
- * Index map for reordering the
- * polynomials.
- */
- std::vector index_map_inverse;
-
- /**
- * Each tensor product polynomial
- * i is a product of
- * one-dimensional polynomials in
- * each space direction. Compute
- * the indices of these
- * one-dimensional polynomials
- * for each space direction,
- * given the index i.
- */
- // fix to avoid compiler warnings about zero
- // length arrays
- void compute_index (const unsigned int i,
- unsigned int (&indices)[(dim>0?dim:1)]) const;
};
/** @} */
+
+/* ---------------- template and inline functions ---------- */
+
#ifndef DOXYGEN
template
+template
inline
-const std::vector &
-TensorProductPolynomialsConst::get_numbering() const
+TensorProductPolynomialsConst::
+TensorProductPolynomialsConst(const std::vector &pols)
+ :
+ TensorProductPolynomials(pols)
{
- return index_map;
+ // append index for renumbering
+ this->index_map.push_back(this->n_tensor_pols);
+ this->index_map_inverse.push_back(this->n_tensor_pols);
}
+
template
inline
-const std::vector &
-TensorProductPolynomialsConst::get_numbering_inverse() const
+unsigned int
+TensorProductPolynomialsConst::n() const
{
- return index_map_inverse;
+ return this->n_tensor_pols+1;
}
-#endif // DOXYGEN
-
-#ifndef DOXYGEN
-/* -------------- declaration of explicit specializations --- */
-
-template <>
-void
-TensorProductPolynomialsConst<1>::compute_index(const unsigned int n,
- unsigned int (&index)[1]) const;
template <>
-void
-TensorProductPolynomialsConst<2>::compute_index(const unsigned int n,
- unsigned int (&index)[2]) const;
-template <>
-void
-TensorProductPolynomialsConst<3>::compute_index(const unsigned int n,
- unsigned int (&index)[3]) const;
-
-
-/* ---------------- template and inline functions ---------- */
-
-template
-template
-TensorProductPolynomialsConst::
-TensorProductPolynomialsConst(const std::vector &pols)
- :
- polynomials (pols.begin(), pols.end()),
- n_tensor_pols(Utilities::fixed_power(pols.size())+1),
- index_map(n_tensor_pols),
- index_map_inverse(n_tensor_pols)
+inline
+unsigned int
+TensorProductPolynomialsConst<0>::n() const
{
- // per default set this index map
- // to identity. This map can be
- // changed by the user through the
- // set_numbering() function
- for (unsigned int i=0; i(&data)+comp);;
+ return *(reinterpret_cast(&data)+comp);
}
/**
diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h
index 765bf03568..dc30e7863a 100644
--- a/deal.II/include/deal.II/dofs/dof_tools.h
+++ b/deal.II/include/deal.II/dofs/dof_tools.h
@@ -954,12 +954,14 @@ namespace DoFTools
*
* Otherwise, if @p face_1 and @p face_2 are not active faces, this
* function loops recursively over the children of @p face_1 and @p face_2.
- *
- * For this to work @p face_1 and @p face_2 must have the same
- * refinement history, i.e. either @p face_1 and @p face_2 must be active
- * faces or must be isotropically refined and have the same number of
- * child faces that recursively obey this rule. (The anisotropic case
- * is not yet implemented.)
+ * If only one of the two faces is active, then we recursively iterate
+ * over the children of the non-active ones and make sure that the
+ * solution function on the refined side equals that on the non-refined
+ * face in much the same way as we enforce hanging node constraints
+ * at places where differently refined cells come together. (However,
+ * unlike hanging nodes, we do not enforce the requirement that there
+ * be only a difference of one refinement level between the two sides
+ * of the domain you would like to be periodic).
*
* This routine only constrains DoFs that are not already constrained.
* If this routine encounters a DoF that already is constrained (for
@@ -967,9 +969,6 @@ namespace DoFTools
* constraint (dofs the entry is constrained to, inhomogeneities) is
* kept and nothing happens.
*
- * Furthermore, no DoFs belonging to (or belonging to any descendant
- * of) @p face_2 get constrained or get marked as being constrained.
- *
* The flags in the @p component_mask (see @ref GlossComponentMask)
* denote which components of the finite element space shall be
* constrained with periodic boundary conditions. If it is left as
@@ -1107,7 +1106,7 @@ namespace DoFTools
* boundary with faces belonging to the second boundary with the help
* of @p orthogonal_equality.
*
- * If this matching is successfull it constrains all DoFs associated
+ * If this matching is successful it constrains all DoFs associated
* with the 'first' boundary to the respective DoFs of the 'second'
* boundary respecting the relative orientation of the two faces.
*
@@ -1117,9 +1116,6 @@ namespace DoFTools
* constraint (dofs the entry is constrained to, inhomogeneities) is
* kept and nothing happens.
*
- * Furthermore, no DoFs belonging to the 'second' boundary get
- * constrained or get marked as being constrained.
- *
* The flags in the last parameter, @p component_mask (see @ref
* GlossComponentMask) denote which components of the finite element space
* shall be constrained with periodic boundary conditions. If it is left
diff --git a/deal.II/include/deal.II/fe/fe.h b/deal.II/include/deal.II/fe/fe.h
index 3059dbcc02..fbf458b9fe 100644
--- a/deal.II/include/deal.II/fe/fe.h
+++ b/deal.II/include/deal.II/fe/fe.h
@@ -343,14 +343,11 @@ class FiniteElement : public Subscriptor,
{
public:
/**
- * Base class for internal data.
- * Adds data for second derivatives to
+ * Base class for internal data. Adds data for second derivatives to
* Mapping::InternalDataBase()
*
- * For information about the
- * general purpose of this class,
- * see the documentation of the
- * base class.
+ * For information about the general purpose of this class, see the
+ * documentation of the base class.
*
* @author Guido Kanschat, 2001
*/
@@ -358,35 +355,24 @@ public:
{
public:
/**
- * Destructor. Needed to
- * avoid memory leaks with
- * difference quotients.
+ * Destructor. Needed to avoid memory leaks with difference quotients.
*/
virtual ~InternalDataBase ();
/**
- * Initialize some pointers
- * used in the computation of
- * second derivatives by
- * finite differencing of
- * gradients.
+ * Initialize some pointers used in the computation of second derivatives
+ * by finite differencing of gradients.
*/
void initialize_2nd (const FiniteElement *element,
const Mapping &mapping,
const Quadrature &quadrature);
/**
- * Storage for FEValues
- * objects needed to
- * approximate second
- * derivatives.
+ * Storage for FEValues objects needed to approximate second derivatives.
*
- * The ordering is p+hx,
- * p+hy, p+hz,
- * p-hx, p-hy,
- * p-hz, where unused
- * entries in lower dimensions
- * are missing.
+ * The ordering is p+hx, p+hy, p+hz, p-hx,
+ * p-hy, p-hz, where unused entries in lower dimensions are
+ * missing.
*/
std::vector*> differences;
};
@@ -400,74 +386,43 @@ public:
const std::vector &nonzero_components);
/**
- * Virtual destructor. Makes sure
- * that pointers to this class
- * are deleted properly.
+ * Virtual destructor. Makes sure that pointers to this class are deleted
+ * properly.
*/
virtual ~FiniteElement ();
/**
- * Return a string that uniquely
- * identifies a finite
- * element. The general
- * convention is that this is the
- * class name, followed by the
- * dimension in angle
- * brackets, and the polynomial
- * degree and whatever else is
- * necessary in parentheses. For
- * example, FE_Q<2>(3) is the
- * value returned for a cubic
- * element in 2d.
+ * Return a string that uniquely identifies a finite element. The general
+ * convention is that this is the class name, followed by the dimension in
+ * angle brackets, and the polynomial degree and whatever else is necessary
+ * in parentheses. For example, FE_Q<2>(3) is the value returned
+ * for a cubic element in 2d.
*
- * Systems of elements have their
- * own naming convention, see the
- * FESystem class.
+ * Systems of elements have their own naming convention, see the FESystem
+ * class.
*/
virtual std::string get_name () const = 0;
/**
- * This operator returns a
- * reference to the present
- * object if the argument given
- * equals to zero. While this
- * does not seem particularly
- * useful, it is helpful in
- * writing code that works with
- * both ::DoFHandler and the hp
- * version hp::DoFHandler, since
- * one can then write code like
- * this:
+ * This operator returns a reference to the present object if the argument
+ * given equals to zero. While this does not seem particularly useful, it is
+ * helpful in writing code that works with both ::DoFHandler and the hp
+ * version hp::DoFHandler, since one can then write code like this:
* @code
* dofs_per_cell
* = dof_handler->get_fe()[cell->active_fe_index()].dofs_per_cell;
* @endcode
*
- * This code doesn't work in both
- * situations without the present
- * operator because
- * DoFHandler::get_fe() returns a
- * finite element, whereas
- * hp::DoFHandler::get_fe()
- * returns a collection of finite
- * elements that doesn't offer a
- * dofs_per_cell
- * member variable: one first has
- * to select which finite element
- * to work on, which is done
- * using the
- * operator[]. Fortunately,
- * cell-@>active_fe_index()
- * also works for non-hp classes
- * and simply returns zero in
- * that case. The present
- * operator[] accepts this zero
- * argument, by returning the
- * finite element with index zero
- * within its collection (that,
- * of course, consists only of
- * the present finite element
- * anyway).
+ * This code doesn't work in both situations without the present operator
+ * because DoFHandler::get_fe() returns a finite element, whereas
+ * hp::DoFHandler::get_fe() returns a collection of finite elements that
+ * doesn't offer a dofs_per_cell member variable: one first has
+ * to select which finite element to work on, which is done using the
+ * operator[]. Fortunately, cell-@>active_fe_index() also works
+ * for non-hp classes and simply returns zero in that case. The present
+ * operator[] accepts this zero argument, by returning the finite element
+ * with index zero within its collection (that, of course, consists only of
+ * the present finite element anyway).
*/
const FiniteElement &operator[] (const unsigned int fe_index) const;
@@ -477,139 +432,75 @@ public:
*/
/**
- * Return the value of the
- * @p ith shape function at the
- * point @p p. @p p is a point
- * on the reference element. If
- * the finite element is
- * vector-valued, then return the
- * value of the only non-zero
- * component of the vector value
- * of this shape function. If the
- * shape function has more than
- * one non-zero component (which
- * we refer to with the term
- * non-primitive), then derived
- * classes implementing this
- * function should throw an
- * exception of type
- * ExcShapeFunctionNotPrimitive. In
- * that case, use the
- * shape_value_component()
- * function.
+ * Return the value of the @p ith shape function at the point @p p. @p p is
+ * a point on the reference element. If the finite element is vector-valued,
+ * then return the value of the only non-zero component of the vector value
+ * of this shape function. If the shape function has more than one non-zero
+ * component (which we refer to with the term non-primitive), then derived
+ * classes implementing this function should throw an exception of type
+ * ExcShapeFunctionNotPrimitive. In that case, use the
+ * shape_value_component() function.
*
- * An
- * ExcUnitShapeValuesDoNotExist
- * is thrown if the shape values
- * of the FiniteElement under
- * consideration depends on the
- * shape of the cell in real
- * space.
+ * An ExcUnitShapeValuesDoNotExist is thrown if the shape values of the
+ * FiniteElement under consideration depends on the shape of the cell in
+ * real space.
*/
virtual double shape_value (const unsigned int i,
const Point &p) const;
/**
- * Just like for shape_value(),
- * but this function will be
- * called when the shape function
- * has more than one non-zero
- * vector component. In that
- * case, this function should
- * return the value of the
- * @p component-th vector
- * component of the @p ith shape
- * function at point @p p.
+ * Just like for shape_value(), but this function will be called when the
+ * shape function has more than one non-zero vector component. In that case,
+ * this function should return the value of the @p component-th vector
+ * component of the @p ith shape function at point @p p.
*/
virtual double shape_value_component (const unsigned int i,
const Point &p,
const unsigned int component) const;
/**
- * Return the gradient of the
- * @p ith shape function at the
- * point @p p. @p p is a point
- * on the reference element, and
- * likewise the gradient is the
- * gradient on the unit cell with
- * respect to unit cell
- * coordinates. If
- * the finite element is
- * vector-valued, then return the
- * value of the only non-zero
- * component of the vector value
- * of this shape function. If the
- * shape function has more than
- * one non-zero component (which
- * we refer to with the term
- * non-primitive), then derived
- * classes implementing this
- * function should throw an
- * exception of type
- * ExcShapeFunctionNotPrimitive. In
- * that case, use the
- * shape_grad_component()
- * function.
+ * Return the gradient of the @p ith shape function at the point @p p. @p p
+ * is a point on the reference element, and likewise the gradient is the
+ * gradient on the unit cell with respect to unit cell coordinates. If the
+ * finite element is vector-valued, then return the value of the only
+ * non-zero component of the vector value of this shape function. If the
+ * shape function has more than one non-zero component (which we refer to
+ * with the term non-primitive), then derived classes implementing this
+ * function should throw an exception of type
+ * ExcShapeFunctionNotPrimitive. In that case, use the
+ * shape_grad_component() function.
*
- * An
- * ExcUnitShapeValuesDoNotExist
- * is thrown if the shape values
- * of the FiniteElement under
- * consideration depends on the
- * shape of the cell in real
- * space.
+ * An ExcUnitShapeValuesDoNotExist is thrown if the shape values of the
+ * FiniteElement under consideration depends on the shape of the cell in
+ * real space.
*/
virtual Tensor<1,dim> shape_grad (const unsigned int i,
const Point &p) const;
/**
- * Just like for shape_grad(),
- * but this function will be
- * called when the shape function
- * has more than one non-zero
- * vector component. In that
- * case, this function should
- * return the gradient of the
- * @p component-th vector
- * component of the @p ith shape
- * function at point @p p.
+ * Just like for shape_grad(), but this function will be called when the
+ * shape function has more than one non-zero vector component. In that case,
+ * this function should return the gradient of the @p component-th vector
+ * component of the @p ith shape function at point @p p.
*/
virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
const Point &p,
const unsigned int component) const;
/**
- * Return the tensor of second
- * derivatives of the @p ith
- * shape function at point @p p
- * on the unit cell. The
- * derivatives are derivatives on
- * the unit cell with respect to
- * unit cell coordinates. If
- * the finite element is
- * vector-valued, then return the
- * value of the only non-zero
- * component of the vector value
- * of this shape function. If the
- * shape function has more than
- * one non-zero component (which
- * we refer to with the term
- * non-primitive), then derived
- * classes implementing this
- * function should throw an
- * exception of type
- * ExcShapeFunctionNotPrimitive. In
- * that case, use the
- * shape_grad_grad_component()
- * function.
+ * Return the tensor of second derivatives of the @p ith shape function at
+ * point @p p on the unit cell. The derivatives are derivatives on the unit
+ * cell with respect to unit cell coordinates. If the finite element is
+ * vector-valued, then return the value of the only non-zero component of
+ * the vector value of this shape function. If the shape function has more
+ * than one non-zero component (which we refer to with the term
+ * non-primitive), then derived classes implementing this function should
+ * throw an exception of type ExcShapeFunctionNotPrimitive. In that case,
+ * use the shape_grad_grad_component() function.
*
- * An
- * ExcUnitShapeValuesDoNotExist
- * is thrown if the shape values
- * of the FiniteElement under
- * consideration depends on the
- * shape of the cell in real
- * space.
+ * An ExcUnitShapeValuesDoNotExist is thrown if the shape values of the
+ * FiniteElement under consideration depends on the shape of the cell in
+ * real space.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
const Point &p) const;
@@ -630,21 +521,14 @@ public:
const Point &p,
const unsigned int component) const;
/**
- * Check for non-zero values on a
- * face in order to optimize out
- * matrix elements.
+ * Check for non-zero values on a face in order to optimize out matrix
+ * elements.
*
- * This function returns
- * @p true, if the shape
- * function @p shape_index has
- * non-zero values on the face
- * @p face_index.
+ * This function returns @p true, if the shape function @p shape_index has
+ * non-zero values on the face @p face_index.
*
- * A default implementation is
- * provided in this basis class
- * which always returns @p
- * true. This is the safe way to
- * go.
+ * A default implementation is provided in this basis class which always
+ * returns @p true. This is the safe way to go.
*/
virtual bool has_support_on_face (const unsigned int shape_index,
const unsigned int face_index) const;
@@ -656,374 +540,214 @@ public:
*/
/**
- * Projection from a fine grid
- * space onto a coarse grid
- * space. If this projection
- * operator is associated with a
- * matrix @p P, then the
- * restriction of this matrix
- * @p P_i to a single child cell
- * is returned here.
- *
- * The matrix @p P is the
- * concatenation or the sum of
- * the cell matrices @p P_i,
- * depending on the
- * #restriction_is_additive_flags. This
- * distinguishes interpolation
- * (concatenation) and projection
- * with respect to scalar
+ * Projection from a fine grid space onto a coarse grid space. If this
+ * projection operator is associated with a matrix @p P, then the
+ * restriction of this matrix @p P_i to a single child cell is returned
+ * here.
+ *
+ * The matrix @p P is the concatenation or the sum of the cell matrices @p
+ * P_i, depending on the #restriction_is_additive_flags. This distinguishes
+ * interpolation (concatenation) and projection with respect to scalar
* products (summation).
*
- * Row and column indices are
- * related to coarse grid and
- * fine grid spaces,
- * respectively, consistent with
- * the definition of the
- * associated operator.
- *
- * If projection matrices are not
- * implemented in the derived
- * finite element class, this
- * function aborts with
- * ExcProjectionVoid. You can
- * check whether this is the case
- * by calling the
- * restriction_is_implemented()
- * or the
- * isotropic_restriction_is_implemented()
- * function.
+ * Row and column indices are related to coarse grid and fine grid spaces,
+ * respectively, consistent with the definition of the associated operator.
+ *
+ * If projection matrices are not implemented in the derived finite element
+ * class, this function aborts with ExcProjectionVoid. You can check whether
+ * this is the case by calling the restriction_is_implemented() or the
+ * isotropic_restriction_is_implemented() function.
*/
- const FullMatrix &
+ virtual const FullMatrix &
get_restriction_matrix (const unsigned int child,
const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const;
/**
* Embedding matrix between grids.
*
- * The identity operator from a
- * coarse grid space into a fine
- * grid space is associated with
- * a matrix @p P. The
- * restriction of this matrix @p P_i to
- * a single child cell is
- * returned here.
- *
- * The matrix @p P is the
- * concatenation, not the sum of
- * the cell matrices
- * @p P_i. That is, if the same
- * non-zero entry j,k exists
- * in in two different child
- * matrices @p P_i, the value
- * should be the same in both
- * matrices and it is copied into
- * the matrix @p P only once.
- *
- * Row and column indices are
- * related to fine grid and
- * coarse grid spaces,
- * respectively, consistent with
- * the definition of the
- * associated operator.
- *
- * These matrices are used by
- * routines assembling the
- * prolongation matrix for
- * multi-level methods. Upon
- * assembling the transfer matrix
- * between cells using this
- * matrix array, zero elements in
- * the prolongation matrix are
- * discarded and will not fill up
- * the transfer matrix.
- *
- * If projection matrices are not
- * implemented in the derived
- * finite element class, this
- * function aborts with
- * ExcEmbeddingVoid. You can
- * check whether this is the case
- * by calling the
- * prolongation_is_implemented()
- * or the
- * isotropic_prolongation_is_implemented()
- * function.
+ * The identity operator from a coarse grid space into a fine grid space is
+ * associated with a matrix @p P. The restriction of this matrix @p P_i to a
+ * single child cell is returned here.
+ *
+ * The matrix @p P is the concatenation, not the sum of the cell matrices @p
+ * P_i. That is, if the same non-zero entry j,k exists in in two
+ * different child matrices @p P_i, the value should be the same in both
+ * matrices and it is copied into the matrix @p P only once.
+ *
+ * Row and column indices are related to fine grid and coarse grid spaces,
+ * respectively, consistent with the definition of the associated operator.
+ *
+ * These matrices are used by routines assembling the prolongation matrix
+ * for multi-level methods. Upon assembling the transfer matrix between
+ * cells using this matrix array, zero elements in the prolongation matrix
+ * are discarded and will not fill up the transfer matrix.
+ *
+ * If projection matrices are not implemented in the derived finite element
+ * class, this function aborts with ExcEmbeddingVoid. You can check whether
+ * this is the case by calling the prolongation_is_implemented() or the
+ * isotropic_prolongation_is_implemented() function.
*/
- const FullMatrix &
+ virtual const FullMatrix &
get_prolongation_matrix (const unsigned int child,
const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const;
/**
- * Return whether this element implements
- * its prolongation matrices. The return
- * value also indicates whether a call to
- * the get_prolongation_matrix()
- * function will generate an error or
- * not.
- *
- * Note, that this function
- * returns true only
- * if the prolongation matrices of
- * the isotropic and all
- * anisotropic refinement cases
- * are implemented. If you are
- * interested in the prolongation
- * matrices for isotropic
- * refinement only, use the
- * isotropic_prolongation_is_implemented
- * function instead.
- *
- * This function is mostly here in order
- * to allow us to write more efficient
- * test programs which we run on all
- * kinds of weird elements, and for which
- * we simply need to exclude certain
- * tests in case something is not
- * implemented. It will in general
- * probably not be a great help in
- * applications, since there is not much
- * one can do if one needs these features
- * and they are not implemented. This
- * function could be used to check
- * whether a call to
- * get_prolongation_matrix() will
- * succeed; however, one then still needs
- * to cope with the lack of information
- * this just expresses.
+ * Return whether this element implements its prolongation matrices. The
+ * return value also indicates whether a call to the
+ * get_prolongation_matrix() function will generate an error or not.
+ *
+ * Note, that this function returns true only if the
+ * prolongation matrices of the isotropic and all anisotropic refinement
+ * cases are implemented. If you are interested in the prolongation matrices
+ * for isotropic refinement only, use the
+ * isotropic_prolongation_is_implemented function instead.
+ *
+ * This function is mostly here in order to allow us to write more efficient
+ * test programs which we run on all kinds of weird elements, and for which
+ * we simply need to exclude certain tests in case something is not
+ * implemented. It will in general probably not be a great help in
+ * applications, since there is not much one can do if one needs these
+ * features and they are not implemented. This function could be used to
+ * check whether a call to get_prolongation_matrix() will succeed;
+ * however, one then still needs to cope with the lack of information this
+ * just expresses.
*/
bool prolongation_is_implemented () const;
/**
- * Return whether this element implements
- * its prolongation matrices for isotropic
- * children. The return value also
- * indicates whether a call to the @p
- * get_prolongation_matrix function will
- * generate an error or not.
- *
- * This function is mostly here in order
- * to allow us to write more efficient
- * test programs which we run on all
- * kinds of weird elements, and for which
- * we simply need to exclude certain
- * tests in case something is not
- * implemented. It will in general
- * probably not be a great help in
- * applications, since there is not much
- * one can do if one needs these features
- * and they are not implemented. This
- * function could be used to check
- * whether a call to
- * get_prolongation_matrix() will
- * succeed; however, one then still needs
- * to cope with the lack of information
- * this just expresses.
+ * Return whether this element implements its prolongation matrices for
+ * isotropic children. The return value also indicates whether a call to the
+ * @p get_prolongation_matrix function will generate an error or not.
+ *
+ * This function is mostly here in order to allow us to write more efficient
+ * test programs which we run on all kinds of weird elements, and for which
+ * we simply need to exclude certain tests in case something is not
+ * implemented. It will in general probably not be a great help in
+ * applications, since there is not much one can do if one needs these
+ * features and they are not implemented. This function could be used to
+ * check whether a call to get_prolongation_matrix() will succeed;
+ * however, one then still needs to cope with the lack of information this
+ * just expresses.
*/
bool isotropic_prolongation_is_implemented () const;
/**
- * Return whether this element implements
- * its restriction matrices. The return
- * value also indicates whether a call to
- * the get_restriction_matrix()
- * function will generate an error or
- * not.
- *
- * Note, that this function
- * returns true only
- * if the restriction matrices of
- * the isotropic and all
- * anisotropic refinement cases
- * are implemented. If you are
- * interested in the restriction
- * matrices for isotropic
- * refinement only, use the
- * isotropic_restriction_is_implemented
- * function instead.
- *
- * This function is mostly here in order
- * to allow us to write more efficient
- * test programs which we run on all
- * kinds of weird elements, and for which
- * we simply need to exclude certain
- * tests in case something is not
- * implemented. It will in general
- * probably not be a great help in
- * applications, since there is not much
- * one can do if one needs these features
- * and they are not implemented. This
- * function could be used to check
- * whether a call to
- * get_restriction_matrix() will
- * succeed; however, one then still needs
- * to cope with the lack of information
- * this just expresses.
+ * Return whether this element implements its restriction matrices. The
+ * return value also indicates whether a call to the
+ * get_restriction_matrix() function will generate an error or not.
+ *
+ * Note, that this function returns true only if the
+ * restriction matrices of the isotropic and all anisotropic refinement
+ * cases are implemented. If you are interested in the restriction matrices
+ * for isotropic refinement only, use the
+ * isotropic_restriction_is_implemented function instead.
+ *
+ * This function is mostly here in order to allow us to write more efficient
+ * test programs which we run on all kinds of weird elements, and for which
+ * we simply need to exclude certain tests in case something is not
+ * implemented. It will in general probably not be a great help in
+ * applications, since there is not much one can do if one needs these
+ * features and they are not implemented. This function could be used to
+ * check whether a call to get_restriction_matrix() will succeed;
+ * however, one then still needs to cope with the lack of information this
+ * just expresses.
*/
bool restriction_is_implemented () const;
/**
- * Return whether this element implements
- * its restriction matrices for isotropic
- * children. The return value also
- * indicates whether a call to the @p
- * get_restriction_matrix function will
- * generate an error or not.
- *
- * This function is mostly here in order
- * to allow us to write more efficient
- * test programs which we run on all
- * kinds of weird elements, and for which
- * we simply need to exclude certain
- * tests in case something is not
- * implemented. It will in general
- * probably not be a great help in
- * applications, since there is not much
- * one can do if one needs these features
- * and they are not implemented. This
- * function could be used to check
- * whether a call to
- * get_restriction_matrix() will
- * succeed; however, one then still needs
- * to cope with the lack of information
- * this just expresses.
+ * Return whether this element implements its restriction matrices for
+ * isotropic children. The return value also indicates whether a call to the
+ * @p get_restriction_matrix function will generate an error or not.
+ *
+ * This function is mostly here in order to allow us to write more efficient
+ * test programs which we run on all kinds of weird elements, and for which
+ * we simply need to exclude certain tests in case something is not
+ * implemented. It will in general probably not be a great help in
+ * applications, since there is not much one can do if one needs these
+ * features and they are not implemented. This function could be used to
+ * check whether a call to get_restriction_matrix() will succeed;
+ * however, one then still needs to cope with the lack of information this
+ * just expresses.
*/
bool isotropic_restriction_is_implemented () const;
/**
- * Access the
- * #restriction_is_additive_flags
- * field. See there for more
+ * Access the #restriction_is_additive_flags field. See there for more
* information on its contents.
*
- * The index must be between zero
- * and the number of shape
- * functions of this element.
+ * The index must be between zero and the number of shape functions of this
+ * element.
*/
bool restriction_is_additive (const unsigned int index) const;
/**
- * Return a readonly reference to
- * the matrix which describes the
- * constraints at the interface
- * between a refined and an
- * unrefined cell.
- *
- * The matrix is obviously empty
- * in only one dimension,
- * since there are no constraints
- * then.
- *
- * Note that some finite elements
- * do not (yet) implement hanging
- * node constraints. If this is
- * the case, then this function
- * will generate an exception,
- * since no useful return value
- * can be generated. If you
- * should have a way to live with
- * this, then you might want to
- * use the
- * constraints_are_implemented()
- * function to check up front
- * whethehr this function will
- * succeed or generate the
- * exception.
+ * Return a read only reference to the matrix which describes the
+ * constraints at the interface between a refined and an unrefined cell.
+ *
+ * The matrix is obviously empty in only one dimension, since there are no
+ * constraints then.
+ *
+ * Note that some finite elements do not (yet) implement hanging node
+ * constraints. If this is the case, then this function will generate an
+ * exception, since no useful return value can be generated. If you should
+ * have a way to live with this, then you might want to use the
+ * constraints_are_implemented() function to check up front whether this
+ * function will succeed or generate the exception.
*/
const FullMatrix &constraints (const dealii::internal::SubfaceCase &subface_case=dealii::internal::SubfaceCase::case_isotropic) const;
/**
- * Return whether this element
- * implements its hanging node
- * constraints. The return value
- * also indicates whether a call
- * to the constraints() function
+ * Return whether this element implements its hanging node constraints. The
+ * return value also indicates whether a call to the constraints() function
* will generate an error or not.
*
- * This function is mostly here
- * in order to allow us to write
- * more efficient test programs
- * which we run on all kinds of
- * weird elements, and for which
- * we simply need to exclude
- * certain tests in case hanging
- * node constraints are not
- * implemented. It will in
- * general probably not be a
- * great help in applications,
- * since there is not much one
- * can do if one needs hanging
- * node constraints and they are
- * not implemented. This function
- * could be used to check whether
- * a call to constraints()
- * will succeed; however, one
- * then still needs to cope with
- * the lack of information this
+ * This function is mostly here in order to allow us to write more efficient
+ * test programs which we run on all kinds of weird elements, and for which
+ * we simply need to exclude certain tests in case hanging node constraints
+ * are not implemented. It will in general probably not be a great help in
+ * applications, since there is not much one can do if one needs hanging
+ * node constraints and they are not implemented. This function could be
+ * used to check whether a call to constraints() will succeed;
+ * however, one then still needs to cope with the lack of information this
* just expresses.
*/
bool constraints_are_implemented (const dealii::internal::SubfaceCase &subface_case=dealii::internal::SubfaceCase::case_isotropic) const;
/**
- * Return whether this element
- * implements its hanging node
- * constraints in the new way,
- * which has to be used to make
- * elements "hp compatible".
- * That means, the element properly
- * implements the
- * get_face_interpolation_matrix
- * and get_subface_interpolation_matrix
- * methods. Therefore the return
- * value also indicates whether a call
- * to the get_face_interpolation_matrix()
- * method and the get_subface_interpolation_matrix()
- * method will generate an error or not.
- *
- * Currently the main purpose of this
- * function is to allow the
- * make_hanging_node_constraints method
- * to decide whether the new procedures,
- * which are supposed to work in the hp
- * framework can be used, or if the old
- * well verified but not hp capable
- * functions should be used. Once the
- * transition to the new scheme for
- * computing the interface constraints is
- * complete, this function will be
- * superfluous and will probably go away.
- *
- * Derived classes should implement this
- * function accordingly. The default
- * assumption is that a finite element
- * does not provide hp capable face
- * interpolation, and the default
- * implementation therefore returns @p
- * false.
+ * Return whether this element implements its hanging node constraints in
+ * the new way, which has to be used to make elements "hp compatible". That
+ * means, the element properly implements the get_face_interpolation_matrix
+ * and get_subface_interpolation_matrix methods. Therefore the return value
+ * also indicates whether a call to the get_face_interpolation_matrix()
+ * method and the get_subface_interpolation_matrix() method will generate an
+ * error or not.
+ *
+ * Currently the main purpose of this function is to allow the
+ * make_hanging_node_constraints method to decide whether the new
+ * procedures, which are supposed to work in the hp framework can be used,
+ * or if the old well verified but not hp capable functions should be used.
+ * Once the transition to the new scheme for computing the interface
+ * constraints is complete, this function will be superfluous and will
+ * probably go away.
+ *
+ * Derived classes should implement this function accordingly. The default
+ * assumption is that a finite element does not provide hp capable face
+ * interpolation, and the default implementation therefore returns @p false.
*/
virtual bool hp_constraints_are_implemented () const;
/**
- * Return the matrix
- * interpolating from the given
- * finite element to the present
- * one. The size of the matrix is
- * then #dofs_per_cell times
+ * Return the matrix interpolating from the given finite element to the
+ * present one. The size of the matrix is then #dofs_per_cell times
* source.#dofs_per_cell.
*
- * Derived elements will have to
- * implement this function. They
- * may only provide interpolation
- * matrices for certain source
- * finite elements, for example
- * those from the same family. If
- * they don't implement
- * interpolation from a given
- * element, then they must throw
- * an exception of type
+ * Derived elements will have to implement this function. They may only
+ * provide interpolation matrices for certain source finite elements, for
+ * example those from the same family. If they don't implement interpolation
+ * from a given element, then they must throw an exception of type
* ExcInterpolationNotImplemented.
*/
virtual void
@@ -1038,24 +762,14 @@ public:
/**
- * Return the matrix
- * interpolating from a face of
- * of one element to the face of
- * the neighboring element.
- * The size of the matrix is
- * then source.#dofs_per_face times
- * this->#dofs_per_face.
+ * Return the matrix interpolating from a face of of one element to the face
+ * of the neighboring element. The size of the matrix is then
+ * source.#dofs_per_face times this->#dofs_per_face.
*
- * Derived elements will have to
- * implement this function. They
- * may only provide interpolation
- * matrices for certain source
- * finite elements, for example
- * those from the same family. If
- * they don't implement
- * interpolation from a given
- * element, then they must throw
- * an exception of type
+ * Derived elements will have to implement this function. They may only
+ * provide interpolation matrices for certain source finite elements, for
+ * example those from the same family. If they don't implement interpolation
+ * from a given element, then they must throw an exception of type
* ExcInterpolationNotImplemented.
*/
virtual void
@@ -1064,24 +778,14 @@ public:
/**
- * Return the matrix
- * interpolating from a face of
- * of one element to the subface of
- * the neighboring element.
- * The size of the matrix is
- * then source.#dofs_per_face times
- * this->#dofs_per_face.
- *
- * Derived elements will have to
- * implement this function. They
- * may only provide interpolation
- * matrices for certain source
- * finite elements, for example
- * those from the same family. If
- * they don't implement
- * interpolation from a given
- * element, then they must throw
- * an exception of type
+ * Return the matrix interpolating from a face of of one element to the
+ * subface of the neighboring element. The size of the matrix is then
+ * source.#dofs_per_face times this->#dofs_per_face.
+ *
+ * Derived elements will have to implement this function. They may only
+ * provide interpolation matrices for certain source finite elements, for
+ * example those from the same family. If they don't implement interpolation
+ * from a given element, then they must throw an exception of type
* ExcInterpolationNotImplemented.
*/
virtual void
@@ -1134,15 +838,11 @@ public:
hp_quad_dof_identities (const FiniteElement &fe_other) const;
/**
- * Return whether this element dominates
- * the one given as argument when they
- * meet at a common face,
- * whether it is the other way around,
- * whether neither dominates, or if
- * either could dominate.
+ * Return whether this element dominates the one given as argument when they
+ * meet at a common face, whether it is the other way around, whether
+ * neither dominates, or if either could dominate.
*
- * For a definition of domination, see
- * FiniteElementBase::Domination and in
+ * For a definition of domination, see FiniteElementBase::Domination and in
* particular the @ref hp_paper "hp paper".
*/
virtual
@@ -1152,21 +852,13 @@ public:
//@}
/**
- * Comparison operator. We also
- * check for equality of the
- * constraint matrix, which is
- * quite an expensive operation.
- * Do therefore use this function
- * with care, if possible only
- * for debugging purposes.
+ * Comparison operator. We also check for equality of the constraint matrix,
+ * which is quite an expensive operation. Do therefore use this function
+ * with care, if possible only for debugging purposes.
*
- * Since this function is not
- * that important, we avoid an
- * implementational question
- * about comparing arrays and do
- * not compare the matrix arrays
- * #restriction and
- * #prolongation.
+ * Since this function is not that important, we avoid an implementational
+ * question about comparing arrays and do not compare the matrix arrays
+ * #restriction and #prolongation.
*/
bool operator == (const FiniteElement &) const;
@@ -1175,116 +867,65 @@ public:
* @{
*/
/**
- * Compute vector component and
- * index of this shape function
- * within the shape functions
- * corresponding to this
- * component from the index of a
- * shape function within this
- * finite element.
+ * Compute vector component and index of this shape function within the
+ * shape functions corresponding to this component from the index of a shape
+ * function within this finite element.
*
- * If the element is scalar, then
- * the component is always zero,
- * and the index within this
- * component is equal to the
- * overall index.
- *
- * If the shape function
- * referenced has more than one
- * non-zero component, then it
- * cannot be associated with one
- * vector component, and an
- * exception of type
- * ExcShapeFunctionNotPrimitive
- * will be raised.
- *
- * Note that if the element is
- * composed of other (base)
- * elements, and a base element
- * has more than one component
- * but all its shape functions
- * are primitive (i.e. are
- * non-zero in only one
- * component), then this mapping
- * contains valid
- * information. However, the
- * index of a shape function of
- * this element within one
- * component (i.e. the second
- * number of the respective entry
- * of this array) does not
- * indicate the index of the
- * respective shape function
- * within the base element (since
- * that has more than one
- * vector-component). For this
- * information, refer to the
- * #system_to_base_table field
- * and the
- * system_to_base_index()
- * function.
+ * If the element is scalar, then the component is always zero, and the
+ * index within this component is equal to the overall index.
+ *
+ * If the shape function referenced has more than one non-zero component,
+ * then it cannot be associated with one vector component, and an exception
+ * of type ExcShapeFunctionNotPrimitive will be raised.
+ *
+ * Note that if the element is composed of other (base) elements, and a base
+ * element has more than one component but all its shape functions are
+ * primitive (i.e. are non-zero in only one component), then this mapping
+ * contains valid information. However, the index of a shape function of
+ * this element within one component (i.e. the second number of the
+ * respective entry of this array) does not indicate the index of the
+ * respective shape function within the base element (since that has more
+ * than one vector-component). For this information, refer to the
+ * #system_to_base_table field and the system_to_base_index() function.
*
- * The use of this function is
- * explained extensively in the
- * step-8 and @ref
- * step_20 "step-20" tutorial
- * programs as well as in the
- * @ref vector_valued module.
+ * The use of this function is explained extensively in the step-8 and @ref
+ * step_20 "step-20" tutorial programs as well as in the @ref vector_valued
+ * module.
*/
std::pair
system_to_component_index (const unsigned int index) const;
/**
- * Compute the shape function for
- * the given vector component and
- * index.
+ * Compute the shape function for the given vector component and index.
*
- * If the element is scalar, then
- * the component must be zero,
- * and the index within this
- * component is equal to the
- * overall index.
+ * If the element is scalar, then the component must be zero, and the index
+ * within this component is equal to the overall index.
*
- * This is the opposite operation
- * from the system_to_component_index()
+ * This is the opposite operation from the system_to_component_index()
* function.
*/
unsigned int component_to_system_index(const unsigned int component,
const unsigned int index) const;
/**
- * Same as
- * system_to_component_index(),
- * but do it for shape functions
- * and their indices on a
- * face. The range of allowed
- * indices is therefore
+ * Same as system_to_component_index(), but do it for shape functions and
+ * their indices on a face. The range of allowed indices is therefore
* 0..#dofs_per_face.
*
- * You will rarely need this
- * function in application
- * programs, since almost all
- * application codes only need to
- * deal with cell indices, not
- * face indices. The function is
- * mainly there for use inside
- * the library.
+ * You will rarely need this function in application programs, since almost
+ * all application codes only need to deal with cell indices, not face
+ * indices. The function is mainly there for use inside the library.
*/
std::pair
face_system_to_component_index (const unsigned int index) const;
/**
- * For faces with non-standard
- * face_orientation in 3D, the dofs on
- * faces (quads) have to be permuted in
- * order to be combined with the correct
- * shape functions. Given a local dof @p
- * index on a quad, return the local index,
- * if the face has non-standard
- * face_orientation, face_flip or
- * face_rotation. In 2D and 1D there is no
- * need for permutation and consequently
- * an exception is thrown.
+ * For faces with non-standard face_orientation in 3D, the dofs on faces
+ * (quads) have to be permuted in order to be combined with the correct
+ * shape functions. Given a local dof @p index on a quad, return the local
+ * index, if the face has non-standard face_orientation, face_flip or
+ * face_rotation. In 2D and 1D there is no need for permutation and
+ * consequently an exception is thrown.
*/
unsigned int adjust_quad_dof_index_for_face_orientation (const unsigned int index,
const bool face_orientation,
@@ -1292,268 +933,166 @@ public:
const bool face_rotation) const;
/**
- * For lines with non-standard
- * line_orientation in 3D, the dofs on
- * lines have to be permuted in order to be
- * combined with the correct shape
- * functions. Given a local dof @p index on
- * a line, return the local index, if the
- * line has non-standard
- * line_orientation. In 2D and 1D there is
- * no need for permutation, so the given
- * index is simply returned.
+ * For lines with non-standard line_orientation in 3D, the dofs on lines
+ * have to be permuted in order to be combined with the correct shape
+ * functions. Given a local dof @p index on a line, return the local index,
+ * if the line has non-standard line_orientation. In 2D and 1D there is no
+ * need for permutation, so the given index is simply returned.
*/
unsigned int adjust_line_dof_index_for_line_orientation (const unsigned int index,
const bool line_orientation) const;
/**
- * Return in which of the vector
- * components of this finite
- * element the @p ith shape
- * function is non-zero. The
- * length of the returned array
- * is equal to the number of
- * vector components of this
- * element.
+ * Return in which of the vector components of this finite element the @p
+ * ith shape function is non-zero. The length of the returned array is equal
+ * to the number of vector components of this element.
+ *
+ * For most finite element spaces, the result of this function will be a
+ * vector with exactly one element being @p true, since for most spaces the
+ * individual vector components are independent. In that case, the component
+ * with the single zero is also the first element of what
+ * system_to_component_index() returns.
*
- * For most finite element
- * spaces, the result of this
- * function will be a vector with
- * exactly one element being
- * @p true, since for most
- * spaces the individual vector
- * components are independent. In
- * that case, the component with
- * the single zero is also the
- * first element of what
- * system_to_component_index()
- * returns.
- *
- * Only for those spaces that couple the
- * components, for example to make a
- * shape function divergence free, will
- * there be more than one @p true entry.
- * Elements for which this is true are
- * called non-primitive (see
- * @ref GlossPrimitive).
+ * Only for those spaces that couple the components, for example to make a
+ * shape function divergence free, will there be more than one @p true
+ * entry. Elements for which this is true are called non-primitive (see
+ * @ref GlossPrimitive).
*/
const ComponentMask &
get_nonzero_components (const unsigned int i) const;
/**
- * Return in how many vector
- * components the @p ith shape
- * function is non-zero. This
- * value equals the number of
- * entries equal to @p true in
- * the result of the
- * get_nonzero_components()
- * function.
+ * Return in how many vector components the @p ith shape function is
+ * non-zero. This value equals the number of entries equal to @p true in the
+ * result of the get_nonzero_components() function.
*
- * For most finite element
- * spaces, the result will be
- * equal to one. It is not equal
- * to one only for those ansatz
- * spaces for which vector-valued
- * shape functions couple the
- * individual components, for
- * example in order to make them
- * divergence-free.
+ * For most finite element spaces, the result will be equal to one. It is
+ * not equal to one only for those ansatz spaces for which vector-valued
+ * shape functions couple the individual components, for example in order to
+ * make them divergence-free.
*/
unsigned int
n_nonzero_components (const unsigned int i) const;
/**
- * Return whether the @p ith
- * shape function is primitive in
- * the sense that the shape
- * function is non-zero in only
- * one vector
- * component. Non-primitive shape
- * functions would then, for
- * example, be those of
- * divergence free ansatz spaces,
- * in which the individual vector
+ * Return whether the @p ith shape function is primitive in the sense that
+ * the shape function is non-zero in only one vector
+ * component. Non-primitive shape functions would then, for example, be
+ * those of divergence free ansatz spaces, in which the individual vector
* components are coupled.
*
- * The result of the function is
- * @p true if and only if the
- * result of
- * n_nonzero_components(i) is
- * equal to one.
+ * The result of the function is @p true if and only if the result of
+ * n_nonzero_components(i) is equal to one.
*/
bool
is_primitive (const unsigned int i) const;
/**
- * Import function that is overloaded
- * by the one above and would otherwise
+ * Import function that is overloaded by the one above and would otherwise
* be hidden.
*/
using FiniteElementData::is_primitive;
/**
- * Number of base elements in a
- * mixed discretization.
- *
- * Note that even for vector
- * valued finite elements, the
- * number of components needs not
- * coincide with the number of
- * base elements, since they may
- * be reused. For example, if you
- * create a FESystem with
- * three identical finite element
- * classes by using the
- * constructor that takes one
- * finite element and a
- * multiplicity, then the number
- * of base elements is still one,
- * although the number of
- * components of the finite
- * element is equal to the
- * multiplicity.
+ * Number of base elements in a mixed discretization.
+ *
+ * Note that even for vector valued finite elements, the number of
+ * components needs not coincide with the number of base elements, since
+ * they may be reused. For example, if you create a FESystem with three
+ * identical finite element classes by using the constructor that takes one
+ * finite element and a multiplicity, then the number of base elements is
+ * still one, although the number of components of the finite element is
+ * equal to the multiplicity.
*/
unsigned int n_base_elements () const;
/**
- * Access to base element
- * objects. If the element is
- * atomic, then
- * base_element(0) is
- * @p this.
+ * Access to base element objects. If the element is atomic, then
+ * base_element(0) is @p this.
*/
virtual
const FiniteElement &
base_element (const unsigned int index) const;
/**
- * This index denotes how often
- * the base element @p index is
- * used in a composed element. If
- * the element is atomic, then
- * the result is always equal to
- * one. See the documentation for
- * the n_base_elements()
- * function for more details.
+ * This index denotes how often the base element @p index is used in a
+ * composed element. If the element is atomic, then the result is always
+ * equal to one. See the documentation for the n_base_elements() function
+ * for more details.
*/
unsigned int
element_multiplicity (const unsigned int index) const;
/**
- * Return for shape function
- * @p index the base element it
- * belongs to, the number of the
- * copy of this base element
- * (which is between zero and the
- * multiplicity of this element),
- * and the index of this shape
- * function within this base
- * element.
+ * Return for shape function @p index the base element it belongs to, the
+ * number of the copy of this base element (which is between zero and the
+ * multiplicity of this element), and the index of this shape function
+ * within this base element.
*
- * If the element is not composed of
- * others, then base and instance
- * are always zero, and the index
- * is equal to the number of the
- * shape function. If the element
- * is composed of single
- * instances of other elements
- * (i.e. all with multiplicity
- * one) all of which are scalar,
- * then base values and dof
- * indices within this element
- * are equal to the
- * #system_to_component_table. It
- * differs only in case the
- * element is composed of other
- * elements and at least one of
- * them is vector-valued itself.
- *
- * This function returns valid
- * values also in the case of
- * vector-valued
- * (i.e. non-primitive) shape
- * functions, in contrast to the
- * system_to_component_index()
- * function.
+ * If the element is not composed of others, then base and instance are
+ * always zero, and the index is equal to the number of the shape
+ * function. If the element is composed of single instances of other
+ * elements (i.e. all with multiplicity one) all of which are scalar, then
+ * base values and dof indices within this element are equal to the
+ * #system_to_component_table. It differs only in case the element is
+ * composed of other elements and at least one of them is vector-valued
+ * itself.
+ *
+ * This function returns valid values also in the case of vector-valued
+ * (i.e. non-primitive) shape functions, in contrast to the
+ * system_to_component_index() function.
*/
std::pair, unsigned int>
system_to_base_index (const unsigned int index) const;
/**
- * Same as
- * system_to_base_index(), but
- * for degrees of freedom located
- * on a face. The range of allowed
- * indices is therefore
- * 0..#dofs_per_face.
+ * Same as system_to_base_index(), but for degrees of freedom located on a
+ * face. The range of allowed indices is therefore 0..#dofs_per_face.
*
- * You will rarely need this
- * function in application
- * programs, since almost all
- * application codes only need to
- * deal with cell indices, not
- * face indices. The function is
- * mainly there for use inside
- * the library.
+ * You will rarely need this function in application programs, since almost
+ * all application codes only need to deal with cell indices, not face
+ * indices. The function is mainly there for use inside the library.
*/
std::pair, unsigned int>
face_system_to_base_index (const unsigned int index) const;
/**
- * Given a base element number,
- * return the first block of a
- * BlockVector it would generate.
+ * Given a base element number, return the first block of a BlockVector it
+ * would generate.
*/
unsigned int first_block_of_base (const unsigned int b) const;
/**
- * For each vector component,
- * return which base
- * element implements this
- * component and which vector
- * component in this base element
- * this is. This information is
- * only of interest for
- * vector-valued finite elements
- * which are composed of several
- * sub-elements. In that case,
- * one may want to obtain
- * information about the element
- * implementing a certain vector
- * component, which can be done
- * using this function and the
- * FESystem::base_element()
- * function.
+ * For each vector component, return which base element implements this
+ * component and which vector component in this base element this is. This
+ * information is only of interest for vector-valued finite elements which
+ * are composed of several sub-elements. In that case, one may want to
+ * obtain information about the element implementing a certain vector
+ * component, which can be done using this function and the
+ * FESystem::base_element() function.
*
- * If this is a scalar finite
- * element, then the return value
- * is always equal to a pair of
- * zeros.
+ * If this is a scalar finite element, then the return value is always equal
+ * to a pair of zeros.
*/
std::pair
component_to_base_index (const unsigned int component) const;
/**
- * Return the base element for
- * this block and the number of
- * the copy of the base element.
+ * Return the base element for this block and the number of the copy of the
+ * base element.
*/
std::pair
block_to_base_index (const unsigned int block) const;
/**
- * The vector block and the index
- * inside the block for this
- * shape function.
+ * The vector block and the index inside the block for this shape function.
*/
std::pair
system_to_block_index (const unsigned int component) const;
/**
- * The vector block for this
- * component.
+ * The vector block for this component.
*/
unsigned int
component_to_block_index (const unsigned int component) const;
@@ -1566,57 +1105,53 @@ public:
*/
/**
- * Return a component mask with as many elements as this
- * object has vector components and of which exactly the
- * one component is true that corresponds to the given
- * argument. See @ref GlossComponentMask "the glossary"
- * for more information.
+ * Return a component mask with as many elements as this object has vector
+ * components and of which exactly the one component is true that
+ * corresponds to the given argument. See @ref GlossComponentMask "the
+ * glossary" for more information.
*
- * @param scalar An object that represents a single scalar
- * vector component of this finite element.
- * @return A component mask that is false in all components
- * except for the one that corresponds to the argument.
+ * @param scalar An object that represents a single scalar vector component
+ * of this finite element.
+ * @return A component mask that is false in all components except for the
+ * one that corresponds to the argument.
*/
ComponentMask
component_mask (const FEValuesExtractors::Scalar &scalar) const;
/**
- * Return a component mask with as many elements as this
- * object has vector components and of which exactly the
- * dim components are true that correspond to the given
- * argument. See @ref GlossComponentMask "the glossary"
- * for more information.
+ * Return a component mask with as many elements as this object has vector
+ * components and of which exactly the dim components are true
+ * that correspond to the given argument. See @ref GlossComponentMask "the
+ * glossary" for more information.
*
- * @param vector An object that represents dim
- * vector components of this finite element.
- * @return A component mask that is false in all components
- * except for the ones that corresponds to the argument.
+ * @param vector An object that represents dim vector components of this
+ * finite element.
+ * @return A component mask that is false in all components except for the
+ * ones that corresponds to the argument.
*/
ComponentMask
component_mask (const FEValuesExtractors::Vector &vector) const;
/**
- * Return a component mask with as many elements as this
- * object has vector components and of which exactly the
- * dim*(dim+1)/2 components are true that
- * correspond to the given argument. See @ref GlossComponentMask "the glossary"
- * for more information.
+ * Return a component mask with as many elements as this object has vector
+ * components and of which exactly the dim*(dim+1)/2 components
+ * are true that correspond to the given argument. See @ref
+ * GlossComponentMask "the glossary" for more information.
*
- * @param sym_tensor An object that represents dim*(dim+1)/2
- * components of this finite element that are jointly to be
- * interpreted as forming a symmetric tensor.
- * @return A component mask that is false in all components
- * except for the ones that corresponds to the argument.
+ * @param sym_tensor An object that represents dim*(dim+1)/2 components of
+ * this finite element that are jointly to be interpreted as forming a
+ * symmetric tensor. @return A component mask that is false in all
+ * components except for the ones that corresponds to the argument.
*/
ComponentMask
component_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;
/**
* Given a block mask (see @ref GlossBlockMask "this glossary entry"),
- * produce a component mask (see @ref GlossComponentMask "this glossary entry")
- * that represents the components that correspond to the blocks selected in
- * the input argument. This is essentially a conversion operator from
- * BlockMask to ComponentMask.
+ * produce a component mask (see @ref GlossComponentMask "this glossary
+ * entry") that represents the components that correspond to the blocks
+ * selected in the input argument. This is essentially a conversion operator
+ * from BlockMask to ComponentMask.
*
* @param block_mask The mask that selects individual blocks of the finite
* element
@@ -1627,84 +1162,80 @@ public:
component_mask (const BlockMask &block_mask) const;
/**
- * Return a block mask with as many elements as this
- * object has blocks and of which exactly the
- * one component is true that corresponds to the given
- * argument. See @ref GlossBlockMask "the glossary"
- * for more information.
+ * Return a block mask with as many elements as this object has blocks and
+ * of which exactly the one component is true that corresponds to the given
+ * argument. See @ref GlossBlockMask "the glossary" for more information.
*
- * @note This function will only succeed if the scalar referenced
- * by the argument encompasses a complete block. In other words,
- * if, for example, you pass an extractor for the single
- * $x$ velocity and this object represents an FE_RaviartThomas
- * object, then the single scalar object you selected is part
- * of a larger block and consequently there is no block mask that
- * would represent it. The function will then produce an exception.
+ * @note This function will only succeed if the scalar referenced by the
+ * argument encompasses a complete block. In other words, if, for example,
+ * you pass an extractor for the single $x$ velocity and this object
+ * represents an FE_RaviartThomas object, then the single scalar object you
+ * selected is part of a larger block and consequently there is no block
+ * mask that would represent it. The function will then produce an
+ * exception.
*
- * @param scalar An object that represents a single scalar
- * vector component of this finite element.
- * @return A component mask that is false in all components
- * except for the one that corresponds to the argument.
+ * @param scalar An object that represents a single scalar vector component
+ * of this finite element.
+ * @return A component mask that is false in all components except for the
+ * one that corresponds to the argument.
*/
BlockMask
block_mask (const FEValuesExtractors::Scalar &scalar) const;
/**
- * Return a component mask with as many elements as this
- * object has vector components and of which exactly the
- * dim components are true that correspond to the given
- * argument. See @ref GlossBlockMask "the glossary"
- * for more information.
+ * Return a component mask with as many elements as this object has vector
+ * components and of which exactly the dim components are true
+ * that correspond to the given argument. See @ref GlossBlockMask "the
+ * glossary" for more information.
*
* @note The same caveat applies as to the version of the function above:
- * The extractor object passed as argument must be so that it corresponds
- * to full blocks and does not split blocks of this element.
+ * The extractor object passed as argument must be so that it corresponds to
+ * full blocks and does not split blocks of this element.
*
- * @param vector An object that represents dim
- * vector components of this finite element.
- * @return A component mask that is false in all components
- * except for the ones that corresponds to the argument.
+ * @param vector An object that represents dim vector components of this
+ * finite element.
+ * @return A component mask that is false in all components except for the
+ * ones that corresponds to the argument.
*/
BlockMask
block_mask (const FEValuesExtractors::Vector &vector) const;
/**
- * Return a component mask with as many elements as this
- * object has vector components and of which exactly the
- * dim*(dim+1)/2 components are true that
- * correspond to the given argument. See @ref GlossBlockMask "the glossary"
- * for more information.
+ * Return a component mask with as many elements as this object has vector
+ * components and of which exactly the dim*(dim+1)/2 components
+ * are true that correspond to the given argument. See @ref GlossBlockMask
+ * "the glossary" for more information.
*
* @note The same caveat applies as to the version of the function above:
- * The extractor object passed as argument must be so that it corresponds
- * to full blocks and does not split blocks of this element.
+ * The extractor object passed as argument must be so that it corresponds to
+ * full blocks and does not split blocks of this element.
*
- * @param sym_tensor An object that represents dim*(dim+1)/2
- * components of this finite element that are jointly to be
- * interpreted as forming a symmetric tensor.
- * @return A component mask that is false in all components
- * except for the ones that corresponds to the argument.
+ * @param sym_tensor An object that represents dim*(dim+1)/2 components of
+ * this finite element that are jointly to be interpreted as forming a
+ * symmetric tensor.
+ * @return A component mask that is false in all components except for the
+ * ones that corresponds to the argument.
*/
BlockMask
block_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;
/**
- * Given a component mask (see @ref GlossComponentMask "this glossary entry"),
- * produce a block mask (see @ref GlossBlockMask "this glossary entry")
- * that represents the blocks that correspond to the components selected in
- * the input argument. This is essentially a conversion operator from
- * ComponentMask to BlockMask.
+ * Given a component mask (see @ref GlossComponentMask "this glossary
+ * entry"), produce a block mask (see @ref GlossBlockMask "this glossary
+ * entry") that represents the blocks that correspond to the components
+ * selected in the input argument. This is essentially a conversion operator
+ * from ComponentMask to BlockMask.
*
- * @note This function will only succeed if the components referenced
- * by the argument encompasses complete blocks. In other words,
- * if, for example, you pass an component mask for the single
- * $x$ velocity and this object represents an FE_RaviartThomas
- * object, then the single component you selected is part
- * of a larger block and consequently there is no block mask that
- * would represent it. The function will then produce an exception.
+ * @note This function will only succeed if the components referenced by the
+ * argument encompasses complete blocks. In other words, if, for example,
+ * you pass an component mask for the single $x$ velocity and this object
+ * represents an FE_RaviartThomas object, then the single component you
+ * selected is part of a larger block and consequently there is no block
+ * mask that would represent it. The function will then produce an
+ * exception.
*
- * @param component_mask The mask that selects individual components of the finite
- * element
+ * @param component_mask The mask that selects individual components of the
+ * finite element
* @return A mask that selects those blocks corresponding to the selected
* blocks of the input argument.
*/
@@ -1719,199 +1250,118 @@ public:
*/
/**
- * Return the support points of
- * the trial functions on the
- * unit cell, if the derived
- * finite element defines some.
- * Finite elements that allow
- * some kind of interpolation
- * operation usually have support
- * points. On the other hand,
- * elements that define their
- * degrees of freedom by, for
- * example, moments on faces, or
- * as derivatives, don't have
- * support points. In that case,
- * the returned field is empty.
- *
- * If the finite element defines
- * support points, then their
- * number equals the number of
- * degrees of freedom of the
- * element. The order of points
- * in the array matches that
- * returned by the
- * cell->get_dof_indices
+ * Return the support points of the trial functions on the unit cell, if the
+ * derived finite element defines some. Finite elements that allow some
+ * kind of interpolation operation usually have support points. On the other
+ * hand, elements that define their degrees of freedom by, for example,
+ * moments on faces, or as derivatives, don't have support points. In that
+ * case, the returned field is empty.
+ *
+ * If the finite element defines support points, then their number equals
+ * the number of degrees of freedom of the element. The order of points in
+ * the array matches that returned by the cell->get_dof_indices
* function.
*
- * See the class documentation
- * for details on support points.
+ * See the class documentation for details on support points.
*/
const std::vector > &
get_unit_support_points () const;
/**
- * Return whether a finite
- * element has defined support
- * points. If the result is true,
- * then a call to the
- * get_unit_support_points()
- * yields a non-empty array.
+ * Return whether a finite element has defined support points. If the result
+ * is true, then a call to the get_unit_support_points() yields a non-empty
+ * array.
*
- * The result may be false if an
- * element is not defined by
- * interpolating shape functions,
- * for example by P-elements on
- * quadrilaterals. It will
- * usually only be true if the
- * element constructs its shape
- * functions by the requirement
- * that they be one at a certain
- * point and zero at all the
- * points associated with the
- * other shape functions.
- *
- * In composed elements (i.e. for
- * the FESystem class, the
- * result will be true if all all
- * the base elements have defined
- * support points.
+ * The result may be false if an element is not defined by interpolating
+ * shape functions, for example by P-elements on quadrilaterals. It will
+ * usually only be true if the element constructs its shape functions by the
+ * requirement that they be one at a certain point and zero at all the
+ * points associated with the other shape functions.
+ *
+ * In composed elements (i.e. for the FESystem class, the result will be
+ * true if all all the base elements have defined support points.
*/
bool has_support_points () const;
/**
- * Return the position of the
- * support point of the
- * @p indexth shape function. If
- * it does not exist, raise an
- * exception.
+ * Return the position of the support point of the @p indexth shape
+ * function. If it does not exist, raise an exception.
*
- * The default implementation
- * simply returns the respective
- * element from the array you get
- * from
- * get_unit_support_points(),
- * but derived elements may
- * overload this function. In
- * particular, note that the
- * FESystem class overloads
- * it so that it can return the
- * support points of individual
- * base elements, if not all the
- * base elements define support
- * points. In this way, you can
- * still ask for certain support
- * points, even if
- * get_unit_support_points()
- * only returns an empty array.
+ * The default implementation simply returns the respective element from the
+ * array you get from get_unit_support_points(), but derived elements may
+ * overload this function. In particular, note that the FESystem class
+ * overloads it so that it can return the support points of individual base
+ * elements, if not all the base elements define support points. In this
+ * way, you can still ask for certain support points, even if
+ * get_unit_support_points() only returns an empty array.
*/
virtual
Point
unit_support_point (const unsigned int index) const;
/**
- * Return the support points of
- * the trial functions on the
- * unit face, if the derived
- * finite element defines some.
- * Finite elements that allow
- * some kind of interpolation
- * operation usually have support
- * points. On the other hand,
- * elements that define their
- * degrees of freedom by, for
- * example, moments on faces, or
- * as derivatives, don't have
- * support points. In that case,
- * the returned field is empty
- *
- * Note that elements that have
- * support points need not
- * necessarily have some on the
- * faces, even if the
- * interpolation points are
- * located physically on a
- * face. For example, the
- * discontinuous elements have
- * interpolation points on the
- * vertices, and for higher
- * degree elements also on the
- * faces, but they are not
- * defined to be on faces since
- * in that case degrees of
- * freedom from both sides of a
- * face (or from all adjacent
- * elements to a vertex) would be
- * identified with each other,
- * which is not what we would
- * like to have). Logically,
- * these degrees of freedom are
- * therefore defined to belong to
- * the cell, rather than the face
- * or vertex. In that case, the
- * returned element would
- * therefore have length zero.
- *
- * If the finite element defines
- * support points, then their
- * number equals the number of
- * degrees of freedom on the face
- * (#dofs_per_face). The order
- * of points in the array matches
- * that returned by the
- * cell->get_dof_indices
- * function.
- *
- * See the class documentation
- * for details on support points.
+ * Return the support points of the trial functions on the unit face, if the
+ * derived finite element defines some. Finite elements that allow some
+ * kind of interpolation operation usually have support points. On the other
+ * hand, elements that define their degrees of freedom by, for example,
+ * moments on faces, or as derivatives, don't have support points. In that
+ * case, the returned field is empty
+ *
+ * Note that elements that have support points need not necessarily have
+ * some on the faces, even if the interpolation points are located
+ * physically on a face. For example, the discontinuous elements have
+ * interpolation points on the vertices, and for higher degree elements also
+ * on the faces, but they are not defined to be on faces since in that case
+ * degrees of freedom from both sides of a face (or from all adjacent
+ * elements to a vertex) would be identified with each other, which is not
+ * what we would like to have). Logically, these degrees of freedom are
+ * therefore defined to belong to the cell, rather than the face or
+ * vertex. In that case, the returned element would therefore have length
+ * zero.
+ *
+ * If the finite element defines support points, then their number equals
+ * the number of degrees of freedom on the face (#dofs_per_face). The order
+ * of points in the array matches that returned by the
+ * cell->get_dof_indices function.
+ *
+ * See the class documentation for details on support points.
*/
const std::vector > &
get_unit_face_support_points () const;
/**
- * Return whether a finite
- * element has defined support
- * points on faces. If the result
- * is true, then a call to the
- * get_unit_face_support_points()
+ * Return whether a finite element has defined support points on faces. If
+ * the result is true, then a call to the get_unit_face_support_points()
* yields a non-empty array.
*
- * For more information, see the
- * documentation for the
- * has_support_points()
+ * For more information, see the documentation for the has_support_points()
* function.
*/
bool has_face_support_points () const;
/**
- * The function corresponding to
- * the unit_support_point()
- * function, but for faces. See
- * there for more information.
+ * The function corresponding to the unit_support_point() function, but for
+ * faces. See there for more information.
*/
virtual
Point
unit_face_support_point (const unsigned int index) const;
/**
- * Return a support point vector
- * for generalized interpolation.
- *
- * See the @ref GlossGeneralizedSupport "glossary entry on generalized support points"
- * for more information.
+ * Return a support point vector for generalized interpolation.
+ *
+ * See the @ref GlossGeneralizedSupport "glossary entry on generalized
+ * support points" for more information.
*/
const std::vector > &
get_generalized_support_points () const;
/**
- * Returns true if the
- * class provides nonempty
- * vectors either from
- * get_unit_support_points() or
- * get_generalized_support_points().
- *
- * See the @ref GlossGeneralizedSupport "glossary entry on generalized support points"
- * for more information.
+ * Returns true if the class provides nonempty vectors either from
+ * get_unit_support_points() or get_generalized_support_points().
+ *
+ * See the @ref GlossGeneralizedSupport "glossary entry on generalized
+ * support points" for more information.
*/
bool has_generalized_support_points () const;
@@ -1922,61 +1372,43 @@ public:
get_generalized_face_support_points () const;
/**
- * Return whether a finite
- * element has defined
- * generalized support
- * points on faces. If the result
- * is true, then a call to the
- * get_generalized_face_support_points
- * yields a non-empty array.
+ * Return whether a finite element has defined generalized support points on
+ * faces. If the result is true, then a call to the
+ * get_generalized_face_support_points yields a non-empty array.
*
- * For more information, see the
- * documentation for the
- * has_support_points()
+ * For more information, see the documentation for the has_support_points()
* function.
*/
bool has_generalized_face_support_points () const;
/**
- * Interpolate a set of scalar
- * values, computed in the
- * generalized support points.
+ * Interpolate a set of scalar values, computed in the generalized support
+ * points.
*
- * @note This function is
- * implemented in
- * FiniteElement for the case
- * that the element has support
- * points. In this case, the
- * resulting coefficients are
- * just the values in the suport
- * points. All other elements
- * must reimplement it.
+ * @note This function is implemented in FiniteElement for the case that the
+ * element has support points. In this case, the resulting coefficients are
+ * just the values in the suport points. All other elements must reimplement
+ * it.
*/
virtual void interpolate(std::vector &local_dofs,
const std::vector &values) const;
/**
- * Interpolate a set of vector
- * values, computed in the
- * generalized support points.
+ * Interpolate a set of vector values, computed in the generalized support
+ * points.
*
- * Since a finite element often
- * only interpolates part of a
- * vector, offset is
- * used to determine the first
- * component of the vector to be
- * interpolated. Maybe consider
- * changing your data structures
- * to use the next function.
+ * Since a finite element often only interpolates part of a vector,
+ * offset is used to determine the first component of the vector to
+ * be interpolated. Maybe consider changing your data structures to use the
+ * next function.
*/
virtual void interpolate(std::vector &local_dofs,
const std::vector > &values,
unsigned int offset = 0) const;
/**
- * Interpolate a set of vector
- * values, computed in the
- * generalized support points.
+ * Interpolate a set of vector values, computed in the generalized support
+ * points.
*/
virtual void interpolate(
std::vector &local_dofs,
@@ -1985,15 +1417,12 @@ public:
//@}
/**
- * Determine an estimate for the
- * memory consumption (in bytes)
- * of this object.
+ * Determine an estimate for the memory consumption (in bytes) of this
+ * object.
*
- * This function is made virtual,
- * since finite element objects
- * are usually accessed through
- * pointers to their base class,
- * rather than the class itself.
+ * This function is made virtual, since finite element objects are usually
+ * accessed through pointers to their base class, rather than the class
+ * itself.
*/
virtual std::size_t memory_consumption () const;
/**
@@ -2023,29 +1452,24 @@ public:
DeclException0 (ExcUnitShapeValuesDoNotExist);
/**
- * Attempt to access support
- * points of a finite element
- * which is not Lagrangian.
+ * Attempt to access support points of a finite element which is not
+ * Lagrangian.
*
* @ingroup Exceptions
*/
DeclException0 (ExcFEHasNoSupportPoints);
/**
- * Attempt to access embedding
- * matrices of a finite element
- * which did not implement these
- * matrices.
+ * Attempt to access embedding matrices of a finite element which did not
+ * implement these matrices.
*
* @ingroup Exceptions
*/
DeclException0 (ExcEmbeddingVoid);
/**
- * Attempt to access restriction
- * matrices of a finite element
- * which did not implement these
- * matrices.
+ * Attempt to access restriction matrices of a finite element which did not
+ * implement these matrices.
*
* Exception
* @ingroup Exceptions
@@ -2053,10 +1477,8 @@ public:
DeclException0 (ExcProjectionVoid);
/**
- * Attempt to access constraint
- * matrices of a finite element
- * which did not implement these
- * matrices.
+ * Attempt to access constraint matrices of a finite element which did not
+ * implement these matrices.
*
* Exception
* @ingroup Exceptions
@@ -2102,235 +1524,140 @@ public:
protected:
/**
- * Reinit the vectors of
- * restriction and prolongation
- * matrices to the right sizes:
- * For every refinement case,
- * except for
- * RefinementCase::no_refinement,
- * and for every child of that
- * refinement case the space of
- * one restriction and
- * prolongation matrix is
- * allocated, see the
- * documentation of the
- * restriction and prolongation
- * vectors for more detail on the
- * actual vector sizes.
- *
- * @param
- * isotropic_restriction_only
- * only the restriction matrices
- * required for isotropic
- * refinement are reinited to the
- * right size.
- * @param
- * isotropic_prolongation_only
- * only the prolongation matrices
- * required for isotropic
- * refinement are reinited to the
- * right size.
+ * Reinit the vectors of restriction and prolongation matrices to the right
+ * sizes: For every refinement case, except for
+ * RefinementCase::no_refinement, and for every child of that refinement
+ * case the space of one restriction and prolongation matrix is allocated,
+ * see the documentation of the restriction and prolongation vectors for
+ * more detail on the actual vector sizes.
+ *
+ * @param isotropic_restriction_only only the restriction matrices required
+ * for isotropic refinement are reinited to the right size.
+ * @param isotropic_prolongation_only only the prolongation matrices
+ * required for isotropic refinement are reinited to the right size.
*/
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false,
const bool isotropic_prolongation_only=false);
/**
- * Vector of projection
- * matrices. See
- * get_restriction_matrix()
- * above. The constructor
- * initializes these matrices to
- * zero dimensions, which can be
- * changed by derived classes
- * implementing them.
- *
- * Note, that
- * restriction[refinement_case-1][child]
- * includes the restriction
- * matrix of child
- * child for the
- * RefinementCase
- * refinement_case. Here,
- * we use
- * refinement_case-1
- * instead of
- * refinement_case
- * as for
- * RefinementCase::no_refinement(=0)
- * there are no restriction
- * matrices available.
+ * Vector of projection matrices. See get_restriction_matrix() above. The
+ * constructor initializes these matrices to zero dimensions, which can be
+ * changed by derived classes implementing them.
+ *
+ * Note, that restriction[refinement_case-1][child] includes
+ * the restriction matrix of child child for the RefinementCase
+ * refinement_case. Here, we use refinement_case-1
+ * instead of refinement_case as for
+ * RefinementCase::no_refinement(=0) there are no restriction matrices
+ * available.
*/
std::vector > > restriction;
/**
- * Vector of embedding
- * matrices. See
- * get_prolongation_matrix()
- * above. The constructor
- * initializes these matrices to
- * zero dimensions, which can be
- * changed by derived classes
- * implementing them.
- *
- * Note, that
- * prolongation[refinement_case-1][child]
- * includes the prolongation
- * matrix of child
- * child for the
- * RefinementCase
- * refinement_case. Here,
- * we use
- * refinement_case-1
- * instead of
- * refinement_case
- * as for
- * RefinementCase::no_refinement(=0)
- * there are no prolongation
- * matrices available.
+ * Vector of embedding matrices. See get_prolongation_matrix()
+ * above. The constructor initializes these matrices to zero dimensions,
+ * which can be changed by derived classes implementing them.
+ *
+ * Note, that prolongation[refinement_case-1][child] includes
+ * the prolongation matrix of child child for the
+ * RefinementCase refinement_case. Here, we use
+ * refinement_case-1 instead of refinement_case as
+ * for RefinementCase::no_refinement(=0) there are no prolongation matrices
+ * available.
*/
std::vector > > prolongation;
/**
- * Specify the constraints which
- * the dofs on the two sides of a
- * cell interface underly if the
- * line connects two cells of
- * which one is refined once.
+ * Specify the constraints which the dofs on the two sides of a cell
+ * interface underly if the line connects two cells of which one is refined
+ * once.
*
- * For further details see the
- * general description of the
- * derived class.
+ * For further details see the general description of the derived class.
*
- * This field is obviously
- * useless in one dimension
- * and has there a zero size.
+ * This field is obviously useless in one dimension and has there a zero
+ * size.
*/
FullMatrix interface_constraints;
/**
- * List of support points on the
- * unit cell, in case the finite
- * element has any. The
- * constructor leaves this field
- * empty, derived classes may
- * write in some contents.
- *
- * Finite elements that allow
- * some kind of interpolation
- * operation usually have support
- * points. On the other hand,
- * elements that define their
- * degrees of freedom by, for
- * example, moments on faces, or
- * as derivatives, don't have
- * support points. In that case,
- * this field remains empty.
+ * List of support points on the unit cell, in case the finite element has
+ * any. The constructor leaves this field empty, derived classes may write
+ * in some contents.
+ *
+ * Finite elements that allow some kind of interpolation operation usually
+ * have support points. On the other hand, elements that define their
+ * degrees of freedom by, for example, moments on faces, or as derivatives,
+ * don't have support points. In that case, this field remains empty.
*/
std::vector > unit_support_points;
/**
- * Same for the faces. See the
- * description of the
- * get_unit_face_support_points()
- * function for a discussion of
- * what contributes a face
- * support point.
+ * Same for the faces. See the description of the
+ * get_unit_face_support_points() function for a discussion of what
+ * contributes a face support point.
*/
std::vector > unit_face_support_points;
/**
- * Support points used for
- * interpolation functions of
- * non-Lagrangian elements.
+ * Support points used for interpolation functions of non-Lagrangian
+ * elements.
*/
std::vector > generalized_support_points;
/**
- * Face support points used for
- * interpolation functions of
- * non-Lagrangian elements.
+ * Face support points used for interpolation functions of non-Lagrangian
+ * elements.
*/
std::vector > generalized_face_support_points;
/**
- * For faces with non-standard
- * face_orientation in 3D, the dofs on
- * faces (quads) have to be permuted in
- * order to be combined with the correct
- * shape functions. Given a local dof @p
- * index on a quad, return the shift in the
- * local index, if the face has
- * non-standard face_orientation,
- * i.e. old_index + shift =
- * new_index. In 2D and 1D there is
- * no need for permutation so the vector is
- * empty. In 3D it has the size of
- * #dofs_per_quad * 8 , where 8 is
- * the number of orientations, a face can
- * be in (all combinations of the three
- * bool flags face_orientation, face_flip
- * and face_rotation).
- *
- * The standard implementation fills this
- * with zeros, i.e. no permuatation at
- * all. Derived finite element classes have
- * to fill this Table with the correct
- * values.
+ * For faces with non-standard face_orientation in 3D, the dofs on faces
+ * (quads) have to be permuted in order to be combined with the correct
+ * shape functions. Given a local dof @p index on a quad, return the shift
+ * in the local index, if the face has non-standard face_orientation,
+ * i.e. old_index + shift = new_index. In 2D and 1D there is no
+ * need for permutation so the vector is empty. In 3D it has the size of
+ * #dofs_per_quad * 8 , where 8 is the number of orientations,
+ * a face can be in (all combinations of the three bool flags
+ * face_orientation, face_flip and face_rotation).
+ *
+ * The standard implementation fills this with zeros, i.e. no permuatation
+ * at all. Derived finite element classes have to fill this Table with the
+ * correct values.
*/
Table<2,int> adjust_quad_dof_index_for_face_orientation_table;
/**
- * For lines with non-standard
- * line_orientation in 3D, the dofs on
- * lines have to be permuted in
- * order to be combined with the correct
- * shape functions. Given a local dof @p
- * index on a line, return the shift in the
- * local index, if the line has
- * non-standard line_orientation,
- * i.e. old_index + shift =
- * new_index. In 2D and 1D there is
- * no need for permutation so the vector is
- * empty. In 3D it has the size of
+ * For lines with non-standard line_orientation in 3D, the dofs on lines
+ * have to be permuted in order to be combined with the correct shape
+ * functions. Given a local dof @p index on a line, return the shift in the
+ * local index, if the line has non-standard line_orientation,
+ * i.e. old_index + shift = new_index. In 2D and 1D there is no
+ * need for permutation so the vector is empty. In 3D it has the size of
* #dofs_per_line.
*
- * The standard implementation fills this
- * with zeros, i.e. no permutation at
- * all. Derived finite element classes have
- * to fill this vector with the correct
- * values.
+ * The standard implementation fills this with zeros, i.e. no permutation at
+ * all. Derived finite element classes have to fill this vector with the
+ * correct values.
*/
std::vector adjust_line_dof_index_for_line_orientation_table;
/**
- * Return the size of interface
- * constraint matrices. Since
- * this is needed in every
- * derived finite element class
- * when initializing their size,
- * it is placed into this
- * function, to avoid having to
- * recompute the
- * dimension-dependent size of
- * these matrices each time.
- *
- * Note that some elements do not
- * implement the interface
- * constraints for certain
- * polynomial degrees. In this
- * case, this function still
- * returns the size these
- * matrices should have when
- * implemented, but the actual
- * matrices are empty.
+ * Return the size of interface constraint matrices. Since this is needed in
+ * every derived finite element class when initializing their size, it is
+ * placed into this function, to avoid having to recompute the
+ * dimension-dependent size of these matrices each time.
+ *
+ * Note that some elements do not implement the interface constraints for
+ * certain polynomial degrees. In this case, this function still returns the
+ * size these matrices should have when implemented, but the actual matrices
+ * are empty.
*/
TableIndices<2>
interface_constraints_size () const;
/**
- * Compute second derivatives by
- * finite differences of
- * gradients.
+ * Compute second derivatives by finite differences of gradients.
*/
void compute_2nd (const Mapping &mapping,
const typename Triangulation::cell_iterator &cell,
@@ -2340,174 +1667,107 @@ protected:
FEValuesData &data) const;
/**
- * Given the pattern of nonzero
- * components for each shape
- * function, compute for each
- * entry how many components are
- * non-zero for each shape
- * function. This function is
- * used in the constructor of
- * this class.
+ * Given the pattern of nonzero components for each shape function, compute
+ * for each entry how many components are non-zero for each shape
+ * function. This function is used in the constructor of this class.
*/
static
std::vector
compute_n_nonzero_components (const std::vector &nonzero_components);
/**
- * Determine the values a finite
- * element should compute on
- * initialization of data for
- * FEValues.
- *
- * Given a set of flags
- * indicating what quantities are
- * requested from a FEValues
- * object, update_once() and
- * update_each() compute which
- * values must really be
- * computed. Then, the
- * fill_*_values functions
- * are called with the result of
- * these.
- *
- * Furthermore, values must be
- * computed either on the unit
- * cell or on the physical
- * cell. For instance, the
- * function values of FE_Q do
- * only depend on the quadrature
- * points on the unit
- * cell. Therefore, this flags
- * will be returned by
- * update_once(). The gradients
- * require computation of the
- * covariant transformation
- * matrix. Therefore,
- * @p update_covariant_transformation
- * and @p update_gradients will
- * be returned by
- * update_each().
- *
- * For an example see the same
- * function in the derived class
- * FE_Q.
+ * Determine the values a finite element should compute on initialization of
+ * data for FEValues.
+ *
+ * Given a set of flags indicating what quantities are requested from a
+ * FEValues object, update_once() and update_each() compute which values
+ * must really be computed. Then, the fill_*_values functions are
+ * called with the result of these.
+ *
+ * Furthermore, values must be computed either on the unit cell or on the
+ * physical cell. For instance, the function values of FE_Q do only depend
+ * on the quadrature points on the unit cell. Therefore, this flags will be
+ * returned by update_once(). The gradients require computation of the
+ * covariant transformation matrix. Therefore, @p
+ * update_covariant_transformation and @p update_gradients will be returned
+ * by update_each().
+ *
+ * For an example see the same function in the derived class FE_Q.
*/
virtual UpdateFlags update_once (const UpdateFlags flags) const = 0;
/**
- * Complementary function for
- * update_once().
+ * Complementary function for update_once().
*
- * While update_once() returns
- * the values to be computed on
- * the unit cell for yielding the
- * required data, this function
- * determines the values that
- * must be recomputed on each
- * cell.
+ * While update_once() returns the values to be computed on the unit cell
+ * for yielding the required data, this function determines the values that
+ * must be recomputed on each cell.
*
- * Refer to update_once() for
- * more details.
+ * Refer to update_once() for more details.
*/
virtual UpdateFlags update_each (const UpdateFlags flags) const = 0;
/**
- * A sort of virtual copy
- * constructor. Some places in
- * the library, for example the
- * constructors of FESystem as
- * well as the hp::FECollection
- * class, need to make copies of
- * finite elements without
- * knowing their exact type. They
- * do so through this function.
+ * A sort of virtual copy constructor. Some places in the library, for
+ * example the constructors of FESystem as well as the hp::FECollection
+ * class, need to make copies of finite elements without knowing their exact
+ * type. They do so through this function.
*/
virtual FiniteElement *clone() const = 0;
private:
/**
- * Store what
- * system_to_component_index()
- * will return.
+ * Store what system_to_component_index() will return.
*/
std::vector< std::pair > system_to_component_table;
/**
- * Map between linear dofs and
- * component dofs on face. This
- * is filled with default values
- * in the constructor, but
- * derived classes will have to
- * overwrite the information if
- * necessary.
+ * Map between linear dofs and component dofs on face. This is filled with
+ * default values in the constructor, but derived classes will have to
+ * overwrite the information if necessary.
*
- * By component, we mean the
- * vector component, not the base
- * element. The information thus
- * makes only sense if a shape
- * function is non-zero in only
+ * By component, we mean the vector component, not the base element. The
+ * information thus makes only sense if a shape function is non-zero in only
* one component.
*/
std::vector< std::pair > face_system_to_component_table;
/**
- * For each shape function, store
- * to which base element and
- * which instance of this base
- * element (in case its
- * multiplicity is greater than
- * one) it belongs, and its index
- * within this base element. If
- * the element is not composed of
- * others, then base and instance
- * are always zero, and the index
- * is equal to the number of the
- * shape function. If the element
- * is composed of single
- * instances of other elements
- * (i.e. all with multiplicity
- * one) all of which are scalar,
- * then base values and dof
- * indices within this element
- * are equal to the
- * #system_to_component_table. It
- * differs only in case the
- * element is composed of other
- * elements and at least one of
- * them is vector-valued itself.
- *
- * This array has valid values
- * also in the case of
- * vector-valued
- * (i.e. non-primitive) shape
- * functions, in contrast to the
+ * For each shape function, store to which base element and which instance
+ * of this base element (in case its multiplicity is greater than one) it
+ * belongs, and its index within this base element. If the element is not
+ * composed of others, then base and instance are always zero, and the index
+ * is equal to the number of the shape function. If the element is composed
+ * of single instances of other elements (i.e. all with multiplicity one)
+ * all of which are scalar, then base values and dof indices within this
+ * element are equal to the #system_to_component_table. It differs only in
+ * case the element is composed of other elements and at least one of them
+ * is vector-valued itself.
+ *
+ * This array has valid values also in the case of vector-valued
+ * (i.e. non-primitive) shape functions, in contrast to the
* #system_to_component_table.
*/
std::vector,unsigned int> >
system_to_base_table;
/**
- * Likewise for the indices on
- * faces.
+ * Likewise for the indices on faces.
*/
std::vector,unsigned int> >
face_system_to_base_table;
/**
- * For each base element, store
- * the number of blocks generated
- * by the base and the first block in a block
- * vector it will generate.
+ * For each base element, store the number of blocks generated by the base
+ * and the first block in a block vector it will generate.
*/
BlockIndices base_to_block_indices;
/**
- * The base element establishing
- * a component.
+ * The base element establishing a component.
*
- * For each component number
- * c, the entries have
- * the following meaning:
+ * For each component number c, the entries have the following
+ * meaning:
*
*
table[c].first.first
*
Number of the base element for c.
@@ -2517,147 +1777,82 @@ private:
*
Multiple of the base element for c.
*
*
- * This variable is set to the
- * correct size by the
- * constructor of this class, but
- * needs to be initialized by
- * derived classes, unless its
- * size is one and the only entry
- * is a zero, which is the case
- * for scalar elements. In that
- * case, the initialization by
- * the base class is sufficient.
+ * This variable is set to the correct size by the constructor of this
+ * class, but needs to be initialized by derived classes, unless its size is
+ * one and the only entry is a zero, which is the case for scalar
+ * elements. In that case, the initialization by the base class is
+ * sufficient.
*/
std::vector, unsigned int> >
component_to_base_table;
/**
- * Projection matrices are
- * concatenated or summed up.
- *
- * This flags decides on how the
- * projection matrices of the
- * children of the same father
- * are put together to one
- * operator. The possible modes
- * are concatenation and
- * summation.
- *
- * If the projection is defined
- * by an interpolation operator,
- * the child matrices are
- * concatenated, i.e. values
- * belonging to the same node
- * functional are identified and
- * enter the interpolated value
- * only once. In this case, the
- * flag must be @p false.
- *
- * For projections with respect
- * to scalar products, the child
- * matrices must be summed up to
- * build the complete matrix. The
- * flag should be @p true.
- *
- * For examples of use of these
- * flags, see the places in the
- * library where it is queried.
- *
- * There is one flag per shape
- * function, indicating whether
- * it belongs to the class of
- * shape functions that are
- * additive in the restriction or
- * not.
- *
- * Note that in previous versions
- * of the library, there was one
- * flag per vector component of
- * the element. This is based on
- * the fact that all the shape
- * functions that belong to the
- * same vector component must
- * necessarily behave in the same
- * way, to make things
- * reasonable. However, the
- * problem is that it is
- * sometimes impossible to query
- * this flag in the vector-valued
- * case: this used to be done
- * with the
- * #system_to_component_index
- * function that returns which
- * vector component a shape
- * function is associated
- * with. The point is that since
- * we now support shape functions
- * that are associated with more
- * than one vector component (for
- * example the shape functions of
- * Raviart-Thomas, or Nedelec
- * elements), that function can
- * no more be used, so it can be
- * difficult to find out which
- * for vector component we would
- * like to query the
+ * Projection matrices are concatenated or summed up.
+ *
+ * This flags decides on how the projection matrices of the children of the
+ * same father are put together to one operator. The possible modes are
+ * concatenation and summation.
+ *
+ * If the projection is defined by an interpolation operator, the child
+ * matrices are concatenated, i.e. values belonging to the same node
+ * functional are identified and enter the interpolated value only once. In
+ * this case, the flag must be @p false.
+ *
+ * For projections with respect to scalar products, the child matrices must
+ * be summed up to build the complete matrix. The flag should be @p true.
+ *
+ * For examples of use of these flags, see the places in the library where
+ * it is queried.
+ *
+ * There is one flag per shape function, indicating whether it belongs to
+ * the class of shape functions that are additive in the restriction or not.
+ *
+ * Note that in previous versions of the library, there was one flag per
+ * vector component of the element. This is based on the fact that all the
+ * shape functions that belong to the same vector component must necessarily
+ * behave in the same way, to make things reasonable. However, the problem
+ * is that it is sometimes impossible to query this flag in the
+ * vector-valued case: this used to be done with the
+ * #system_to_component_index function that returns which vector component a
+ * shape function is associated with. The point is that since we now support
+ * shape functions that are associated with more than one vector component
+ * (for example the shape functions of Raviart-Thomas, or Nedelec elements),
+ * that function can no more be used, so it can be difficult to find out
+ * which for vector component we would like to query the
* restriction-is-additive flags.
*/
const std::vector restriction_is_additive_flags;
/**
- * For each shape function, give
- * a vector of bools (with size
- * equal to the number of vector
- * components which this finite
- * element has) indicating in
- * which component each of these
- * shape functions is non-zero.
+ * For each shape function, give a vector of bools (with size equal to the
+ * number of vector components which this finite element has) indicating in
+ * which component each of these shape functions is non-zero.
*
- * For primitive elements, there
- * is only one non-zero
- * component.
+ * For primitive elements, there is only one non-zero component.
*/
const std::vector nonzero_components;
/**
- * This array holds how many
- * values in the respective entry
- * of the #nonzero_components
- * element are non-zero. The
- * array is thus a short-cut to
- * allow faster access to this
- * information than if we had to
- * count the non-zero entries
- * upon each request for this
- * information. The field is
- * initialized in the constructor
- * of this class.
+ * This array holds how many values in the respective entry of the
+ * #nonzero_components element are non-zero. The array is thus a short-cut
+ * to allow faster access to this information than if we had to count the
+ * non-zero entries upon each request for this information. The field is
+ * initialized in the constructor of this class.
*/
const std::vector n_nonzero_components_table;
/**
- * Second derivatives of shapes
- * functions are not computed
- * analytically, but by finite
- * differences of the
- * gradients. This static
- * variable denotes the step
- * length to be used for
- * that. It's value is set to
- * 1e-6.
+ * Second derivatives of shapes functions are not computed analytically, but
+ * by finite differences of the gradients. This static variable denotes the
+ * step length to be used for that. It's value is set to 1e-6.
*/
static const double fd_step_length;
/**
- * Prepare internal data
- * structures and fill in values
- * independent of the
- * cell. Returns a pointer to an
- * object of which the caller of
- * this function then has to
- * assume ownership (which
- * includes destruction when it
- * is no more needed).
+ * Prepare internal data structures and fill in values independent of the
+ * cell. Returns a pointer to an object of which the caller of this function
+ * then has to assume ownership (which includes destruction when it is no
+ * more needed).
*/
virtual typename Mapping::InternalDataBase *
get_data (const UpdateFlags flags,
@@ -2665,16 +1860,10 @@ private:
const Quadrature &quadrature) const = 0;
/**
- * Prepare internal data
- * structure for transformation
- * of faces and fill in values
- * independent of the
- * cell. Returns a pointer to an
- * object of which the caller of
- * this function then has to
- * assume ownership (which
- * includes destruction when it
- * is no more needed).
+ * Prepare internal data structure for transformation of faces and fill in
+ * values independent of the cell. Returns a pointer to an object of which
+ * the caller of this function then has to assume ownership (which includes
+ * destruction when it is no more needed).
*/
virtual typename Mapping::InternalDataBase *
get_face_data (const UpdateFlags flags,
@@ -2682,16 +1871,10 @@ private:
const Quadrature &quadrature) const;
/**
- * Prepare internal data
- * structure for transformation
- * of children of faces and fill
- * in values independent of the
- * cell. Returns a pointer to an
- * object of which the caller of
- * this function then has to
- * assume ownership (which
- * includes destruction when it
- * is no more needed).
+ * Prepare internal data structure for transformation of children of faces
+ * and fill in values independent of the cell. Returns a pointer to an
+ * object of which the caller of this function then has to assume ownership
+ * (which includes destruction when it is no more needed).
*/
virtual typename Mapping::InternalDataBase *
get_subface_data (const UpdateFlags flags,
@@ -2699,15 +1882,11 @@ private:
const Quadrature &quadrature) const;
/**
- * Fill the fields of
- * FEValues. This function
- * performs all the operations
- * needed to compute the data of an
- * FEValues object.
+ * Fill the fields of FEValues. This function performs all the operations
+ * needed to compute the data of an FEValues object.
*
- * The same function in
- * @p mapping must have been
- * called for the same cell first!
+ * The same function in @p mapping must have been called for the same cell
+ * first!
*/
virtual void
fill_fe_values (const Mapping &mapping,
@@ -2719,15 +1898,11 @@ private:
CellSimilarity::Similarity &cell_similarity) const = 0;
/**
- * Fill the fields of
- * FEFaceValues. This function
- * performs all the operations
- * needed to compute the data of an
- * FEFaceValues object.
+ * Fill the fields of FEFaceValues. This function performs all the
+ * operations needed to compute the data of an FEFaceValues object.
*
- * The same function in
- * @p mapping must have been
- * called for the same cell first!
+ * The same function in @p mapping must have been called for the same cell
+ * first!
*/
virtual void
fill_fe_face_values (const Mapping &mapping,
@@ -2739,15 +1914,11 @@ private:
FEValuesData &data) const = 0;
/**
- * Fill the fields of
- * FESubfaceValues. This function
- * performs all the operations
- * needed to compute the data of an
- * FESubfaceValues object.
+ * Fill the fields of FESubfaceValues. This function performs all the
+ * operations needed to compute the data of an FESubfaceValues object.
*
- * The same function in
- * @p mapping must have been
- * called for the same cell first!
+ * The same function in @p mapping must have been called for the same cell
+ * first!
*/
virtual void
fill_fe_subface_values (const Mapping &mapping,
diff --git a/deal.II/include/deal.II/fe/fe_face.h b/deal.II/include/deal.II/fe/fe_face.h
index 8747b6d4f9..4fdddd2e07 100644
--- a/deal.II/include/deal.II/fe/fe_face.h
+++ b/deal.II/include/deal.II/fe/fe_face.h
@@ -47,49 +47,35 @@ class FE_FaceQ : public FE_PolyFace, dim, spaced
{
public:
/**
- * Constructor for tensor product
- * polynomials of degree
- * p. The shape
- * functions created using this
- * constructor correspond to
- * Legendre polynomials in each
- * coordinate direction.
+ * Constructor for tensor product polynomials of degree p. The
+ * shape functions created using this constructor correspond to Lagrange
+ * polynomials in each coordinate direction.
*/
FE_FaceQ(unsigned int p);
virtual FiniteElement *clone() const;
/**
- * Return a string that uniquely
- * identifies a finite
- * element. This class returns
- * FE_FaceQ(degree), with
- * dim and degree
- * replaced by appropriate
- * values.
+ * Return a string that uniquely identifies a finite element. This class
+ * returns FE_FaceQ(degree), with dim and
+ * degree replaced by appropriate values.
*/
virtual std::string get_name () const;
/**
* Check for non-zero values on a face.
*
- * This function returns
- * @p true, if the shape
- * function @p shape_index has
- * non-zero values on the face
- * @p face_index.
+ * This function returns @p true, if the shape function @p shape_index has
+ * non-zero values on the face @p face_index.
*
- * Implementation of the
- * interface in
- * FiniteElement
+ * Implementation of the interface in FiniteElement
*/
virtual bool has_support_on_face (const unsigned int shape_index,
const unsigned int face_index) const;
private:
/**
- * Return vector with dofs per
- * vertex, line, quad, hex.
+ * Return vector with dofs per vertex, line, quad, hex.
*/
static std::vector get_dpo_vector (const unsigned int deg);
};
diff --git a/deal.II/include/deal.II/fe/fe_poly.h b/deal.II/include/deal.II/fe/fe_poly.h
index ba8ccd7297..1ff0992513 100644
--- a/deal.II/include/deal.II/fe/fe_poly.h
+++ b/deal.II/include/deal.II/fe/fe_poly.h
@@ -74,122 +74,85 @@ public:
const std::vector &nonzero_components);
/**
- * Return the polynomial degree
- * of this finite element,
- * i.e. the value passed to the
- * constructor.
+ * Return the polynomial degree of this finite element, i.e. the value
+ * passed to the constructor.
*/
unsigned int get_degree () const;
/**
- * Return the numbering of the underlying
- * polynomial space compared to
- * lexicographic ordering of the basis
- * functions. Returns
+ * Return the numbering of the underlying polynomial space compared to
+ * lexicographic ordering of the basis functions. Returns
* POLY::get_numbering().
*/
std::vector get_poly_space_numbering() const;
/**
- * Return the inverse numbering of the
- * underlying polynomial space. Returns
+ * Return the inverse numbering of the underlying polynomial space. Returns
* POLY::get_numbering_inverse().
*/
std::vector get_poly_space_numbering_inverse() const;
/**
- * Return the value of the
- * ith shape function at
- * the point p. See the
- * FiniteElement base class
- * for more information about the
- * semantics of this function.
+ * Return the value of the ith shape function at the point
+ * p. See the FiniteElement base class for more information about
+ * the semantics of this function.
*/
virtual double shape_value (const unsigned int i,
const Point &p) const;
/**
- * Return the value of the
- * componentth vector
- * component of the ith
- * shape function at the point
- * p. See the
- * FiniteElement base class
- * for more information about the
- * semantics of this function.
+ * Return the value of the componentth vector component of the
+ * ith shape function at the point p. See the
+ * FiniteElement base class for more information about the semantics of this
+ * function.
*
- * Since this element is scalar,
- * the returned value is the same
- * as if the function without the
- * _component suffix
- * were called, provided that the
- * specified component is zero.
+ * Since this element is scalar, the returned value is the same as if the
+ * function without the _component suffix were called, provided
+ * that the specified component is zero.
*/
virtual double shape_value_component (const unsigned int i,
const Point &p,
const unsigned int component) const;
/**
- * Return the gradient of the
- * ith shape function at
- * the point p. See the
- * FiniteElement base class
- * for more information about the
- * semantics of this function.
+ * Return the gradient of the ith shape function at the point
+ * p. See the FiniteElement base class for more information about
+ * the semantics of this function.
*/
virtual Tensor<1,dim> shape_grad (const unsigned int i,
const Point &p) const;
/**
- * Return the gradient of the
- * componentth vector
- * component of the ith
- * shape function at the point
- * p. See the
- * FiniteElement base class
- * for more information about the
- * semantics of this function.
+ * Return the gradient of the componentth vector component of the
+ * ith shape function at the point p. See the
+ * FiniteElement base class for more information about the semantics of this
+ * function.
*
- * Since this element is scalar,
- * the returned value is the same
- * as if the function without the
- * _component suffix
- * were called, provided that the
- * specified component is zero.
+ * Since this element is scalar, the returned value is the same as if the
+ * function without the _component suffix were called, provided
+ * that the specified component is zero.
*/
virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
const Point &p,
const unsigned int component) const;
/**
- * Return the tensor of second
- * derivatives of the
- * ith shape function at
- * point p on the unit
- * cell. See the
- * FiniteElement base class
- * for more information about the
- * semantics of this function.
+ * Return the tensor of second derivatives of the ith shape
+ * function at point p on the unit cell. See the FiniteElement base
+ * class for more information about the semantics of this function.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
const Point &p) const;
/**
- * Return the second derivative
- * of the componentth
- * vector component of the
- * ith shape function at
- * the point p. See the
- * FiniteElement base class
- * for more information about the
- * semantics of this function.
+ * Return the second derivative of the componentth vector component
+ * of the ith shape function at the point p. See the
+ * FiniteElement base class for more information about the semantics of this
+ * function.
*
- * Since this element is scalar,
- * the returned value is the same
- * as if the function without the
- * _component suffix
- * were called, provided that the
- * specified component is zero.
+ * Since this element is scalar, the returned value is the same as if the
+ * function without the _component suffix were called, provided
+ * that the specified component is zero.
*/
virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
const Point &p,
@@ -232,84 +195,46 @@ protected:
/**
- * Determine the values that need
- * to be computed on the unit
- * cell to be able to compute all
- * values required by
- * flags.
+ * Determine the values that need to be computed on the unit cell to be able
+ * to compute all values required by flags.
*
- * For the purpuse of this
- * function, refer to the
- * documentation in
+ * For the purpuse of this function, refer to the documentation in
* FiniteElement.
*
- * This class assumes that shape
- * functions of this
- * FiniteElement do not
- * depend on the actual shape of
- * the cells in real
- * space. Therefore, the effect
- * in this element is as follows:
- * if update_values is
- * set in flags, copy it
- * to the result. All other flags
- * of the result are cleared,
- * since everything else must be
- * computed for each cell.
+ * This class assumes that shape functions of this FiniteElement do
+ * not depend on the actual shape of the cells in real
+ * space. Therefore, the effect in this element is as follows: if
+ * update_values is set in flags, copy it to the
+ * result. All other flags of the result are cleared, since everything else
+ * must be computed for each cell.
*/
virtual UpdateFlags update_once (const UpdateFlags flags) const;
/**
- * Determine the values that need
- * to be computed on every cell
- * to be able to compute all
- * values required by
- * flags.
+ * Determine the values that need to be computed on every cell to be able to
+ * compute all values required by flags.
*
- * For the purpuse of this
- * function, refer to the
- * documentation in
+ * For the purpuse of this function, refer to the documentation in
* FiniteElement.
*
- * This class assumes that shape
- * functions of this
- * FiniteElement do not
- * depend on the actual shape of
- * the cells in real
- * space.
+ * This class assumes that shape functions of this FiniteElement do
+ * not depend on the actual shape of the cells in real space.
*
- * The effect in this element is
- * as follows:
+ * The effect in this element is as follows:
*
-
- *
if
- * update_gradients is
- * set, the result will contain
- * update_gradients and
- * update_covariant_transformation.
- * The latter is required to
- * transform the gradient on the
- * unit cell to the real
- * cell. Remark, that the action
- * required by
- * update_covariant_transformation
- * is actually performed by the
- * Mapping object used in
- * conjunction with this finite
- * element.
if
- * update_hessians
- * is set, the result will
- * contain
- * update_hessians
- * and
- * update_covariant_transformation.
- * The rationale is the same as
- * above and no higher
- * derivatives of the
- * transformation are required,
- * since we use difference
- * quotients for the actual
- * computation.
+ *
+ *
if update_gradients is set, the result will contain
+ * update_gradients and update_covariant_transformation.
+ * The latter is required to transform the gradient on the unit cell to the
+ * real cell. Remark, that the action required by
+ * update_covariant_transformation is actually performed by the
+ * Mapping object used in conjunction with this finite element.
+ *
+ *
if update_hessians is set, the result will contain
+ * update_hessians and update_covariant_transformation.
+ * The rationale is the same as above and no higher derivatives of the
+ * transformation are required, since we use difference quotients for the
+ * actual computation.
*
*
*/
@@ -319,61 +244,37 @@ protected:
/**
* Fields of cell-independent data.
*
- * For information about the
- * general purpose of this class,
- * see the documentation of the
- * base class.
+ * For information about the general purpose of this class, see the
+ * documentation of the base class.
*/
class InternalData : public FiniteElement::InternalDataBase
{
public:
/**
- * Array with shape function
- * values in quadrature
- * points. There is one
- * row for each shape
- * function, containing
- * values for each quadrature
- * point.
+ * Array with shape function values in quadrature points. There is one row
+ * for each shape function, containing values for each quadrature point.
*
- * In this array, we store
- * the values of the shape
- * function in the quadrature
- * points on the unit
- * cell. Since these values
- * do not change under
- * transformation to the real
- * cell, we only need to copy
- * them over when visiting a
- * concrete cell.
+ * In this array, we store the values of the shape function in the
+ * quadrature points on the unit cell. Since these values do not change
+ * under transformation to the real cell, we only need to copy them over
+ * when visiting a concrete cell.
*/
std::vector > shape_values;
/**
- * Array with shape function
- * gradients in quadrature
- * points. There is one
- * row for each shape
- * function, containing
- * values for each quadrature
+ * Array with shape function gradients in quadrature points. There is one
+ * row for each shape function, containing values for each quadrature
* point.
*
- * We store the gradients in
- * the quadrature points on
- * the unit cell. We then
- * only have to apply the
- * transformation (which is a
- * matrix-vector
- * multiplication) when
- * visiting an actual cell.
+ * We store the gradients in the quadrature points on the unit cell. We
+ * then only have to apply the transformation (which is a matrix-vector
+ * multiplication) when visiting an actual cell.
*/
std::vector > > shape_gradients;
};
/**
- * The polynomial space. Its type
- * is given by the template
- * parameter POLY.
+ * The polynomial space. Its type is given by the template parameter POLY.
*/
POLY poly_space;
};
diff --git a/deal.II/include/deal.II/fe/fe_poly.templates.h b/deal.II/include/deal.II/fe/fe_poly.templates.h
index 99f49d7a02..f5ba25953e 100644
--- a/deal.II/include/deal.II/fe/fe_poly.templates.h
+++ b/deal.II/include/deal.II/fe/fe_poly.templates.h
@@ -460,18 +460,18 @@ namespace internal
return std::vector();
}
- template
+ template
inline
std::vector
- get_poly_space_numbering (const TensorProductPolynomials &poly)
+ get_poly_space_numbering (const TensorProductPolynomials &poly)
{
return poly.get_numbering();
}
- template
+ template
inline
std::vector
- get_poly_space_numbering_inverse (const TensorProductPolynomials &poly)
+ get_poly_space_numbering_inverse (const TensorProductPolynomials &poly)
{
return poly.get_numbering_inverse();
}
diff --git a/deal.II/include/deal.II/fe/fe_q.h b/deal.II/include/deal.II/fe/fe_q.h
index e17e1e4e07..294daad459 100644
--- a/deal.II/include/deal.II/fe/fe_q.h
+++ b/deal.II/include/deal.II/fe/fe_q.h
@@ -14,7 +14,7 @@
#include
#include
-#include
+#include
DEAL_II_NAMESPACE_OPEN
@@ -220,386 +220,55 @@ DEAL_II_NAMESPACE_OPEN
* @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann, 2001, 2004, 2005; Oliver Kayser-Herold, 2004; Katharina Kormann, 2008; Martin Kronbichler, 2008
*/
template
-class FE_Q : public FE_Poly,dim,spacedim>
+class FE_Q : public FE_Q_Base,dim,spacedim>
{
public:
/**
- * Constructor for tensor product
- * polynomials of degree @p p.
+ * Constructor for tensor product polynomials of degree @p p.
*/
FE_Q (const unsigned int p);
/**
- * Constructor for tensor product
- * polynomials with support
- * points @p points based on a
- * one-dimensional quadrature
- * formula. The degree of the
- * finite element is
- * points.size()-1.
- * Note that the first point has
- * to be 0 and the last one 1. If
- * FE_Q(QGaussLobatto<1>(fe_degree+1))
- * is specified, so-called
- * Gauss-Lobatto elements are
- * obtained which can give a
- * diagonal mass matrix if
- * combined with Gauss-Lobatto
- * quadrature on the same
- * points. Their use is shown in
- * step-48.
+ * Constructor for tensor product polynomials with support points @p points
+ * based on a one-dimensional quadrature formula. The degree of the finite
+ * element is points.size()-1. Note that the first point has to be
+ * 0 and the last one 1. If
+ * FE_Q(QGaussLobatto<1>(fe_degree+1)) is specified, so-called
+ * Gauss-Lobatto elements are obtained which can give a diagonal mass matrix
+ * if combined with Gauss-Lobatto quadrature on the same points. Their use
+ * is shown in step-48.
*/
FE_Q (const Quadrature<1> &points);
/**
- * Return a string that uniquely
- * identifies a finite
- * element. This class returns
- * FE_Q(degree), with
- * @p dim and @p degree
- * replaced by appropriate
- * values.
+ * Constructs a FE_Q_isoQ1 element. That element shares large parts of code
+ * with FE_Q so most of the construction work is done in this routine,
+ * whereas the public constructor is in the class FE_Q_isoQ1.
*/
- virtual std::string get_name () const;
-
- /**
- * Return the matrix
- * interpolating from the given
- * finite element to the present
- * one. The size of the matrix is
- * then @p dofs_per_cell times
- * source.dofs_per_cell.
- *
- * These matrices are only
- * available if the source
- * element is also a @p FE_Q
- * element. Otherwise, an
- * exception of type
- * FiniteElement::ExcInterpolationNotImplemented
- * is thrown.
- */
- virtual void
- get_interpolation_matrix (const FiniteElement &source,
- FullMatrix &matrix) const;
-
-
- /**
- * Return the matrix
- * interpolating from a face of
- * of one element to the face of
- * the neighboring element.
- * The size of the matrix is
- * then source.dofs_per_face times
- * this->dofs_per_face.
- *
- * Derived elements will have to
- * implement this function. They
- * may only provide interpolation
- * matrices for certain source
- * finite elements, for example
- * those from the same family. If
- * they don't implement
- * interpolation from a given
- * element, then they must throw
- * an exception of type
- * FiniteElement::ExcInterpolationNotImplemented.
- */
- virtual void
- get_face_interpolation_matrix (const FiniteElement &source,
- FullMatrix &matrix) const;
-
- /**
- * Return the matrix
- * interpolating from a face of
- * of one element to the face of
- * the neighboring element.
- * The size of the matrix is
- * then source.dofs_per_face times
- * this->dofs_per_face.
- *
- * Derived elements will have to
- * implement this function. They
- * may only provide interpolation
- * matrices for certain source
- * finite elements, for example
- * those from the same family. If
- * they don't implement
- * interpolation from a given
- * element, then they must throw
- * an exception of type
- * FiniteElement::ExcInterpolationNotImplemented.
- */
- virtual void
- get_subface_interpolation_matrix (const FiniteElement &source,
- const unsigned int subface,
- FullMatrix &matrix) const;
-
- /**
- * Check for non-zero values on a face.
- *
- * This function returns
- * @p true, if the shape
- * function @p shape_index has
- * non-zero values on the face
- * @p face_index.
- *
- * Implementation of the
- * interface in
- * FiniteElement
- */
- virtual bool has_support_on_face (const unsigned int shape_index,
- const unsigned int face_index) const;
-
- /**
- * @name Functions to support hp
- * @{
- */
-
- /**
- * Return whether this element
- * implements its hanging node
- * constraints in the new way,
- * which has to be used to make
- * elements "hp compatible".
- *
- * For the FE_Q class the result is
- * always true (independent of the degree
- * of the element), as it implements the
- * complete set of functions necessary
- * for hp capability.
- */
- virtual bool hp_constraints_are_implemented () const;
-
- /**
- * If, on a vertex, several finite elements are active, the hp code
- * first assigns the degrees of freedom of each of these FEs
- * different global indices. It then calls this function to find out
- * which of them should get identical values, and consequently can
- * receive the same global DoF index. This function therefore
- * returns a list of identities between DoFs of the present finite
- * element object with the DoFs of @p fe_other, which is a reference
- * to a finite element object representing one of the other finite
- * elements active on this particular vertex. The function computes
- * which of the degrees of freedom of the two finite element objects
- * are equivalent, both numbered between zero and the corresponding
- * value of dofs_per_vertex of the two finite elements. The first
- * index of each pair denotes one of the vertex dofs of the present
- * element, whereas the second is the corresponding index of the
- * other finite element.
- */
- virtual
- std::vector >
- hp_vertex_dof_identities (const FiniteElement &fe_other) const;
+ FE_Q(const unsigned int subdivisions_per_dimension,
+ const unsigned int base_degree);
/**
- * Same as hp_vertex_dof_indices(), except that the function treats
- * degrees of freedom on lines.
+ * Return a string that uniquely identifies a finite element. This class
+ * returns FE_Q(degree), with @p dim and @p degree replaced by
+ * appropriate values.
*/
- virtual
- std::vector >
- hp_line_dof_identities (const FiniteElement &fe_other) const;
-
- /**
- * Same as hp_vertex_dof_indices(), except that the function treats
- * degrees of freedom on quads.
- */
- virtual
- std::vector >
- hp_quad_dof_identities (const FiniteElement