Hermite polynomials #14657

ivweber · 2023-01-10T12:11:29Z

I have grouped together the modified files into different sets. These are the files for generating Hermite polynomials.

bangerth

Yes, very nice! I've left a bunch of comments, but this is close -- nicely done, many thanks!

bangerth · 2023-01-12T22:01:39Z

include/deal.II/base/polynomials_hermite.h

+#include <deal.II/lac/full_matrix.h>
+
+#include <algorithm>
+#include <exception>
+#include <memory>
+#include <vector>


I don't think you actually need any of these header files.

One of the member functions returns an std::vector, but you're right I think the others might be from an older implementation.

bangerth · 2023-01-12T22:02:48Z

include/deal.II/base/polynomials_hermite.h

+   * either a non-zero shape value or derivative at $y=0$ and $y=1$ 
+   * on the reference interval $y \in [0,1]$.


This is an interesting convention to use $y$ to denote the single variable the polynomial depends on. Would you mind changing it to $x$ here and everywhere below?

bangerth · 2023-01-13T20:04:16Z

include/deal.II/base/polynomials_hermite.h

+   * $j-(r+1)$ (or value for $j=r+1$) at $x=1$. The basis is rescaled
+   * so that a function corresponding to a non-zero $j^{th}$
+   * derivative has derivative value $j! 4^{j}$ at the corresponding
+   * node.


Is it correct that the value of the 0th and (r+1)st basis function is one at the end points of the interval then? Maybe state that as well, just to be extra clear.

Yes, this is correct. I will specify this directly.

bangerth · 2023-01-13T20:05:17Z

include/deal.II/base/polynomials_hermite.h

+   * derivative has derivative value $j! 4^{j}$ at the corresponding
+   * node.
+   */
+  class HermitePoly : public Polynomial<double>


All of the other classes in this direction have names of the form PolynomialsX. So maybe call this one PolynomialsHermite?

bangerth · 2023-01-13T20:16:16Z

include/deal.II/base/polynomials_hermite.h

+  {
+  public:
+    /**
+     * Constructor for an individual Hermite polynomial. We write f_{j}


Suggested change

* Constructor for an individual Hermite polynomial. We write f_{j}

* Constructor for an individual Hermite polynomial. We write $f_{j}$

source/base/polynomials_hermite.cc

bangerth · 2023-01-13T20:26:11Z

source/base/polynomials_hermite.cc

+    , degree(2 * regularity + 1)
+    , regularity(regularity)
+    , side_index(index % (regularity + 1))
+    , side((index > regularity) ? 1 : 0)


This may be easier to read as follows, given the modulus operation in the line above:

Suggested change

, side((index > regularity) ? 1 : 0)

, side((index >= regularity+1) ? 1 : 0)

bangerth · 2023-01-13T20:27:10Z

source/base/polynomials_hermite.cc

+    , regularity(regularity)
+    , side_index(index % (regularity + 1))
+    , side((index > regularity) ? 1 : 0)
+  {}


Would you mind adding an assertion like this:

Suggested change

{}

{

AssertIndexRange(index, 2*(regularity+1));

}

bangerth · 2023-01-13T20:27:45Z

source/base/polynomials_hermite.cc

+    const unsigned int              sz = 2 * regularity + 2;
+
+    for (unsigned int i = 0; i < sz; ++i)
+      polys.push_back(HermitePoly(regularity, i));


Suggested change

polys.push_back(HermitePoly(regularity, i));

polys.emplace_back(HermitePoly(regularity, i));

source/base/polynomials_hermite.cc

peterrum

@ivweber Looks good to me once you fixed the few typos and you have fixed the merge conflict!

@drwells @bangerth Could you also give a final look at the PR so that Ivy can proceed with the follow-up PRs.

peterrum · 2023-05-13T06:44:13Z

.gitignore

@@ -16,6 +16,7 @@ CMakeLists.txt.user
 /compile_commands.json
 /.clang_complete
 /contrib/utilities/programs
+/include/build/


Suggested change

/include/build/

peterrum · 2023-05-13T06:45:43Z

include/deal.II/base/polynomials_hermite.h

+  /**
+   * This class implements Hermite interpolation polynomials (see
+   * @cite CiarletRiavart1972interpolation) enforcing the maximum
+   * posible level of regularity $r$ in the FEM basis given a


Suggested change

* posible level of regularity $r$ in the FEM basis given a

* possible level of regularity $r$ in the FEM basis given a

peterrum · 2023-05-13T06:47:19Z

include/deal.II/base/polynomials_hermite.h

+    /**
+     * This stores whether the shape function corresponds to a non-zero
+     * value or derivative at $x=0$ on the reference interval
+     * ($\mathtt{side} =0$), or at $x=1$ ($\mathtt{side} =1$).


Suggested change

* ($\mathtt{side} =0$), or at $x=1$ ($\mathtt{side} =1$).

* ($\mathtt{side} =0$) or at $x=1$ ($\mathtt{side} =1$).

ivweber · 2023-05-15T10:47:26Z

@peterrum I've just corrected the typos and changed .gitignore back.

drwells

This looks good but we could improve things by picking a clearer name.

include/deal.II/base/polynomials_hermite.h

peterrum · 2023-05-15T21:08:01Z

@ivweber May I ask you to resolve the conflict in doc/doxygen/references.bib and to rebase the branch on master. It seems that some tests are not running which might be because the changes to the tests have been integrated in the meantime. You could optionally also squash the commits.

ivweber · 2023-05-16T07:53:28Z

@peterrum @drwells I have just rebased the pull request and corrected the references file.

ivweber · 2023-05-16T08:50:11Z

I've just reverted back to a clean version from before the rebase

kronbichler

Thank you for tackling this. To me, both the implementation and the naming look good, as the first step to implement a finite element based on these polynomials.

Just to be sure: I assume it would be difficult/impossible to express these polynomials as the product of the factors (x-x0)*(x-x0)*(x-x_1)*(x-x_1) for some x_0 and x_1? For the values it is easy to see what to do, but not for the derivatives. The reason I'm asking is that for Lagrange polynomials we observed that expressing the polynomial with conventional coefficients in front of 1, x, x^2, etc. are unstable in terms of roundoff behavior, and one needs to form with roots. I'm wondering if the Hermite polynomials have similar behavior. On the other hand, it really starts to matter for degrees beyond 6, so we might be fine. If you're wondering and have not looked into it, the evaluation stability is the reason why the Polynomial class has these members:

dealii/include/deal.II/base/polynomial.h

Lines 304 to 314 in cb6ee41

    
               /** 
        
                * 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$. 
        
                */ 
        
               std::vector<number> lagrange_support_points;

And for evaluation

dealii/include/deal.II/base/polynomial.h

Lines 819 to 839 in cb6ee41

    
           if (in_lagrange_product_form == false) 
        
             { 
        
               Assert(coefficients.size() > 0, ExcEmptyObject()); 
        
               // Horner scheme 
        
               const unsigned int m     = coefficients.size(); 
        
               number             value = coefficients.back(); 
        
               for (int k = m - 2; k >= 0; --k) 
        
                 value = value * x + coefficients[k]; 
        
               return value; 
        
             } 
        
           else 
        
             { 
        
               // direct evaluation of Lagrange polynomial 
        
               const unsigned int m     = lagrange_support_points.size(); 
        
               number             value = 1.; 
        
               for (unsigned int j = 0; j < m; ++j) 
        
                 value *= x - lagrange_support_points[j]; 
        
               value *= lagrange_weight; 
        
               return value; 
        
             }

kronbichler · 2023-05-16T12:48:21Z

include/deal.II/base/polynomials_hermite.h

+#ifndef dealii_hermite_polynomials
+#define dealii_hermite_polynomials


Can you please rename:

Suggested change

#ifndef dealii_hermite_polynomials

#define dealii_hermite_polynomials

#ifndef dealii_polynomials_hermite_h

#define dealii_polynomials_hermite_h

peterrum · 2023-05-22T08:24:43Z

@ivweber I think we are very close to merging!

However, some of the tests are not picked up. Only 4 tests are picked up in contrast to the 10, which should be the normal case (see #15165). I think the might be related to the fact that the branch is not rebased on master (not 100% sure). May I ask you to rebase the code. If it helps I can do that for you and I force push it to this branch.

ivweber · 2023-05-22T08:35:19Z

@kronbichler As far as I know, for most Hermite interpolation polynomials it's difficult to find a fully factorised form of the polynomial because part of the expression is a truncated Taylor series at x=0 for (1-x)^(-r-1), where r is the level of regularity being enforced. Some of the polynomials can be easily written like this, but it's only for the whole basis at lower polynomial degrees.

kronbichler

Looks good, apart from a minor comment regarding appearance. Thank you for the explanation, this is what I suspected.

kronbichler · 2023-05-22T09:11:58Z

doc/doxygen/references.bib

@@ -2159,4 +2168,4 @@ @article{zalesak1979fully
  pages={335--362},
  year={1979},
  publisher={Elsevier}
-}
+}


Can you please fix the missing line break at the end of the file, such that this diff here disappears?

ivweber · 2023-05-22T09:12:50Z

@peterrum I tried git rebase dealii/master on a copy of the branch and removed the duplicate entries in references.bib that this produced, but I can still only see 4 tests. Is it possible that github-actions doesn't have access to this PR for some reason?

kronbichler · 2023-05-22T09:16:34Z

/rebuild

peterrum · 2023-05-22T09:18:35Z

I tried git rebase dealii/master on a copy of the branch and removed the duplicate entries in references.bib that this produced, but I can still only see 4 tests. Is it possible that github-actions doesn't have access to this PR for some reason?

I must admit I am bit clueless what the issue is. Maybe you need to create a new PR!? This would be really annoying...

ivweber · 2023-05-25T08:52:12Z

@drwells Is there a way to get GitHub Actions to start tests on this PR? I looked at the missing tests mentioned above, and they were all ones directly linked to GitHub. Also for some reason the indent test isn't listed as a GitHub Action here, while it is on other PRs.

drwells · 2023-05-25T14:29:54Z

I don't know - it should work. I'll spend some time digging into this problem.

drwells · 2023-05-25T17:50:39Z

I think we ran into an interesting corner-case where the CI needs to be approved for first time contributors, but the first PR (#13239) wasn't actually merged or had CI enabled. I need someone with higher repository permissions to handle this (or possibly upgrade my access).

ivweber · 2023-05-30T09:27:02Z

@drwells Would it be useful to remove the previous pull request, or would be likely to create further issues?

tamiko · 2023-06-05T23:58:11Z

@ivweber I have pushed a rebase of your changes to your branch (on github) in order to get the tests going.

peterrum · 2023-06-06T19:53:46Z

@ivweber Finally, the test have run. There are few indentation issues. Good to go once you fixed these 👍

tamiko · 2023-06-09T17:37:47Z

@ivweber I will rebase this pull request again to get rid of your merge commit. This means that the history on github and on your local branch will diverge.

Edit: I have squashed everything into a single commit to clean up the history and making it easier to review the changes.

bangerth · 2023-06-11T00:24:37Z

Nice job, @ivweber ! It took a little while to get this merged, but welcome to the deal.II authors!

drwells self-assigned this Jan 12, 2023

bangerth requested changes Jan 13, 2023

View reviewed changes

ivweber force-pushed the pull_request_polynomials branch from fbe78e3 to bd2904f Compare May 9, 2023 08:01

peterrum approved these changes May 13, 2023

View reviewed changes

peterrum added Finite Element ready to test labels May 13, 2023

ivweber force-pushed the pull_request_polynomials branch from bd2904f to 715f9a7 Compare May 15, 2023 10:43

peterrum approved these changes May 15, 2023

View reviewed changes

drwells approved these changes May 15, 2023

View reviewed changes

include/deal.II/base/polynomials_hermite.h Show resolved Hide resolved

ivweber force-pushed the pull_request_polynomials branch 4 times, most recently from 8c2c02a to 876b8e8 Compare May 16, 2023 07:45

ivweber force-pushed the pull_request_polynomials branch from 876b8e8 to 8c2c02a Compare May 16, 2023 08:49

kronbichler approved these changes May 16, 2023

View reviewed changes

ivweber force-pushed the pull_request_polynomials branch from 40cbebe to ed3f570 Compare May 22, 2023 08:14

ivweber force-pushed the pull_request_polynomials branch 2 times, most recently from 6a15bb6 to 7b58a85 Compare May 22, 2023 09:08

kronbichler approved these changes May 22, 2023

View reviewed changes

ivweber force-pushed the pull_request_polynomials branch from 7b58a85 to 8ddb800 Compare May 22, 2023 09:14

ivweber requested a review from bangerth May 25, 2023 08:47

bangerth approved these changes May 27, 2023

View reviewed changes

tamiko force-pushed the pull_request_polynomials branch from 8ddb800 to 0715177 Compare June 5, 2023 23:57

base: implement Hermite polynomials

16c8693

tamiko force-pushed the pull_request_polynomials branch from 73ce577 to 16c8693 Compare June 9, 2023 17:56

tamiko merged commit 77cf9aa into dealii:master Jun 10, 2023
14 checks passed

		* either a non-zero shape value or derivative at $y=0$ and $y=1$
		* on the reference interval $y \in [0,1]$.

	* Constructor for an individual Hermite polynomial. We write f_{j}
	* Constructor for an individual Hermite polynomial. We write $f_{j}$

	, side((index > regularity) ? 1 : 0)
	, side((index >= regularity+1) ? 1 : 0)

-  {}
+  {
+    AssertIndexRange(index, 2*(regularity+1));
+  }

	polys.push_back(HermitePoly(regularity, i));
	polys.emplace_back(HermitePoly(regularity, i));

	* posible level of regularity $r$ in the FEM basis given a
	* possible level of regularity $r$ in the FEM basis given a

	* ($\mathtt{side} =0$), or at $x=1$ ($\mathtt{side} =1$).
	* ($\mathtt{side} =0$) or at $x=1$ ($\mathtt{side} =1$).

	/**
	* 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$.
	*/
	std::vector<number> lagrange_support_points;

	if (in_lagrange_product_form == false)
	{
	Assert(coefficients.size() > 0, ExcEmptyObject());

	// Horner scheme
	const unsigned int m = coefficients.size();
	number value = coefficients.back();
	for (int k = m - 2; k >= 0; --k)
	value = value * x + coefficients[k];
	return value;
	}
	else
	{
	// direct evaluation of Lagrange polynomial
	const unsigned int m = lagrange_support_points.size();
	number value = 1.;
	for (unsigned int j = 0; j < m; ++j)
	value *= x - lagrange_support_points[j];
	value *= lagrange_weight;
	return value;
	}

		#ifndef dealii_hermite_polynomials
		#define dealii_hermite_polynomials

Hermite polynomials #14657

Hermite polynomials #14657

Conversation

ivweber commented Jan 10, 2023

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ivweber commented May 15, 2023

drwells left a comment

Choose a reason for hiding this comment

peterrum commented May 15, 2023

ivweber commented May 16, 2023

ivweber commented May 16, 2023

kronbichler left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum commented May 22, 2023

ivweber commented May 22, 2023

kronbichler left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ivweber commented May 22, 2023

kronbichler commented May 22, 2023

peterrum commented May 22, 2023

ivweber commented May 25, 2023

drwells commented May 25, 2023

drwells commented May 25, 2023

ivweber commented May 30, 2023

tamiko commented Jun 5, 2023

peterrum commented Jun 6, 2023

tamiko commented Jun 9, 2023 • edited

bangerth commented Jun 11, 2023

tamiko commented Jun 9, 2023 •

edited