Implicit integration for elastodynamics using Newton's method

General Discussion
1

Following many examples provided in the FreeFem documentation (Non-linear elasticity, Non-linear static problems, Newton Method for the Steady Navier-Stokes equation, Incompressible Navier Stokes, Steady Navier Stokes, etc.), I went on to write a script which tries to perform an implicit integration in time of a beam subject to a mass at its tip using Newton’s method.
More precisely, the numerical example is the same as the one presented in Section 3.7.1 of this PhD manuscript, using a midpoint scheme (see Section 3.4.4 of the same manuscript). The system of equations is defined in Section 3.2.2, equation (3.5).

Issues.
The integration scheme cannot reproduce their numerical results. My vertical displacement is not reaching values larger than 10^{-14} meters. Since I am very new to FreeFem, I think the major source of mistakes lies between the chair and the computer. Even if the code is somewhat complicated, maybe just a glance at it from experts here would suffice to know what I am doing wrong.
I include the theory for this problem just after, how to perform Newton’s method in this case as well as the relevant computations (differentials for the most part), and then the whole code at the bottom as well as some personal comments on it. In particular, I am not sure whether I did the correct construction to apply surface forces only on a part of the mesh. I also hope Newton’s update are done correctly. Any feedback, even minor, is appreciated.

Theory and context.
Denoting by \varphi the displacement, \dot{\varphi} the velocity and p the pressure (seen here as a Lagrange multiplier for the quasi-incompressibility constraint), the midpoint scheme reads

\left\lbrace \begin{aligned} &\int_{\Omega}\rho\dot{\varphi}_{n+1}\cdot v=\int_{\Omega}\rho\dot{\varphi}_n\cdot v-\Delta t_n\int_{\Omega}\Pi_{n+1/2}:\nabla v+\Delta t_n\int_{\Omega}\frac{f_n+f_{n+1}}{2}\cdot v,\\ &\int_{\Omega}\varphi_{n+1}\cdot w=\int_{\Omega}\varphi_n\cdot w+\Delta t_n\int_{\Omega}\frac{\dot{\varphi}_n+\dot{\varphi}_{n+1}}{2}\cdot w,\\ &\int_{\Omega}q\left(D_{n+1/2}-1+\varepsilon\frac{p_n+p_{n+1}}{2}\right)=0. \end{aligned} \right.

where v,w,q are test functions, and f\equiv (0,-f) is a constant force which models the force applied to the beam by a ball on a string (for instance). Here,

\left\lbrace \begin{aligned} \Pi_{n+1/2}&=\partial_F\hat{W}\left(\frac{F_n+F_{n+1}}{2}\right)-\frac{p_n+p_{n+1}}{2}{\rm cof}\left(\frac{F_n+F_{n+1}}{2}\right),\\ D_{n+1/2}&={\rm det}\left(\frac{F_n+F_{n+1}}{2}\right). \end{aligned} \right.

where {\rm cof}(A) denotes the cofactor matrix of a matrix A, and F_n\equiv \nabla\varphi_n. The elastic energy is given by the Mooney-Rivlin law:

\hat{W}(F)=W(C)=c_1({\rm Tr}C-3)+c_2({\rm Tr\, cof} C-3),

where C=F^TF. After some computations (and if no mistakes from my part are present), this simplifies to

\hat{W}(F)=(c_1+c_2)(\left\lVert F\right\rVert^2-3),

where \left\lVert\cdot\right\rVert is the 2-norm (Frobenius norm).

Numerical parameters c_1,c_2,\varepsilon are the same as in the manuscript example.

Newton’s method.
We fix \Delta t_n=\Delta t and f_n=f constants. For a given configuration (\dot{\varphi}^{n},\varphi^{n},p^{n}) at time n\Delta t, we find the new configuration at iteration n+1 using Newton’s method: define F_n(\dot{\varphi},\varphi,p) as

\begin{equation*} F_n(\dot{\varphi},\varphi,p)= \left\lbrace \begin{aligned} &\int_{\Omega}\rho\dot{\varphi}\cdot v-\rho\dot{\varphi}^{n}\cdot v+\Delta t(c_1+c_2)(\nabla\varphi+\nabla\varphi^{n}):\nabla v\\ &-\int_{\Omega}\Delta t\frac{p+p^{n}}{2}{\rm cof}\left(\frac{\nabla\varphi+\nabla\varphi^{n}}{2}\right):\nabla v-\Delta t f\cdot v\\ &+\int_{\Omega}\varphi\cdot w-\varphi^{n}\cdot w-\Delta t\frac{\dot{\varphi}+\dot{\varphi}^{n}}{2}\cdot w\\ &+\int_{\Omega}q\left(\det\left(\frac{\nabla\varphi+\nabla\varphi^{n}}{2}\right)-1+\varepsilon\frac{p+p^{n}}{2}\right) \end{aligned} \right\rbrace. \end{equation*}

We want to find (\dot{\varphi}^{n+1},\varphi^{n+1},p^{n+1}) such that F_n(\dot{\varphi}^{n+1},\varphi^{n+1},p^{n+1})=0. To compute the differential of F_n, we need to expand the two terms involving the cofactor matrix and the determinant. As we are set in dimension 2, the cofactor function is in fact linear, so that

\begin{equation*} {\rm cof}\left(\frac{\nabla\varphi+\nabla\varphi^{n}}{2}\right)=\frac{1}{2}\left({\rm cof}(\nabla\varphi)+{\rm cof}(\nabla\varphi^{n})\right). \end{equation*}

The differential of (\varphi,p)\mapsto\frac{p+p^{n}}{2}{\rm cof}\left(\frac{\nabla\varphi+\nabla\varphi^{n}}{2}\right) at a point (\varphi,p) therefore acts on (\delta\varphi,\delta p) as

\begin{equation*} \frac{1}{4}\left(p^{n}{\rm cof}(\nabla\delta\varphi)+\delta p{\rm cof}(\nabla\varphi^{n})+p{\rm cof}(\nabla\delta\varphi)+\delta p{\rm cof}(\nabla\varphi)\right). \end{equation*}

We now expand the determinant:

\begin{align*} \det\left( \frac{\nabla\varphi+\nabla\varphi^{n}}{2} \right) &=\frac{1}{4} \left[ (\partial_x\varphi_x+\partial_x\varphi_x^{n})(\partial_y\varphi_y+\partial_y\varphi_y^{n})-(\partial_x\varphi_y+\partial_x\varphi_y^{n})(\partial_y\varphi_x+\partial_y\varphi_x^{n}) \right]\\ &= \frac{1}{4}\left[ \partial_x\varphi_x\partial_y\varphi_y +\partial_x\varphi_x\partial_y\varphi_y^{n} +\partial_x\varphi_x^{n}\partial_y\varphi_y +\partial_x\varphi_x^{n}\partial_y\varphi_y^{n} \right]\\ &\quad -\frac{1}{4}\left[ \partial_x\varphi_y\partial_y\varphi_x +\partial_x\varphi_y\partial_y\varphi_x^{n} +\partial_x\varphi_y^{n}\partial_y\varphi_x +\partial_x\varphi_y^{n}\partial_y\varphi_x^{n} \right]. \end{align*}

Therefore, the differential of \varphi\mapsto {\rm det}\left(\frac{\nabla\varphi+\nabla\varphi^{n}}{2}\right) at a point \varphi acts on \delta\varphi as

\begin{align*} &\frac{1}{4}\left[ \partial_x\varphi_x\partial_y\delta\varphi_y+\partial_x\delta\varphi_x\partial_y\varphi_y +\partial_x\delta\varphi_x\partial_y\varphi_y^{n} +\partial_x\varphi_x^{n}\partial_y\delta\varphi_y \right]\\ &\quad-\frac{1}{4}\left[ \partial_x\varphi_y\partial_y\delta\varphi_x+\partial_x\delta\varphi_y\partial_y\varphi_x +\partial_x\delta\varphi_y\partial_y\varphi_x^{n} +\partial_x\varphi_y^{n}\partial_y\delta\varphi_x \right]. \end{align*}

Eventually, the differential of F_n at a point (\dot{\varphi},\varphi,p) acts on (\delta\dot{\varphi},\delta\varphi,\delta p) as

\begin{equation*} DF_n(\dot{\varphi},\varphi,p)(\delta\dot{\varphi},\delta\varphi,\delta p) =\left\lbrace \begin{aligned} &\int_{\Omega}\rho\delta\dot{\varphi}\cdot v+\Delta t(c_1+c_2)\nabla\delta\varphi:\nabla v\\ &-\int_{\Omega}\frac{\Delta t}{4}\left(p^{n}{\rm cof}(\nabla\delta\varphi)+\delta p{\rm cof}(\nabla\varphi^{n})+p{\rm cof}(\nabla\delta\varphi)+\delta p{\rm cof}(\nabla\varphi)\right):\nabla v\\ &+\int_{\Omega}\delta\varphi\cdot w-\frac{\Delta t}{2}\delta\dot{\varphi}\cdot w\\ &+\int_{\Omega}\frac{q}{4}\left[ \partial_x\varphi_x\partial_y\delta\varphi_y+\partial_x\delta\varphi_x\partial_y\varphi_y +\partial_x\delta\varphi_x\partial_y\varphi_y^{n} +\partial_x\varphi_x^{n}\partial_y\delta\varphi_y \right]\\ &-\int_{\Omega}\frac{q}{4}\left[ \partial_x\varphi_y\partial_y\delta\varphi_x+\partial_x\delta\varphi_y\partial_y\varphi_x +\partial_x\delta\varphi_y\partial_y\varphi_x^{n} +\partial_x\varphi_y^{n}\partial_y\delta\varphi_x \right]\\ &+\int_{\Omega}\frac{q}{2}\varepsilon\delta p \end{aligned} \right\rbrace. \end{equation*}

Newton’s method then reads: set an initial condition, e.g. (\dot{\varphi}^{n+1,0},\varphi^{n+1,0},p^{n+1,0})=(\dot{\varphi}^{n},\varphi^{n},p^{n}), and for 0\leqslant i< N_{\rm Newton}, solve

\begin{equation*} DF_n(\dot{\varphi}^{n+1,i},\varphi^{n+1,i},p^{n+1,i})(\delta\dot{\varphi}^{n+1,i},\delta\varphi^{n+1,i},\delta p^{n+1,i})=F_n(\dot{\varphi}^{n+1,i},\varphi^{n+1,i},p^{n+1,i}) \end{equation*}

and update (\dot{\varphi}^{n+1,i+1},\varphi^{n+1,i+1},p^{n+1,i+1})=(\dot{\varphi}^{n+1,i},\varphi^{n+1,i},p^{n+1,i})-(\delta\dot{\varphi}^{n+1,i},\delta\varphi^{n+1,i},\delta p^{n+1,i}); and break if convergence is reached (by checking if the residue, for instance the L^{2}/H^{1} norm, is sufficiently small).

Code.

// Parameters

real length=1;
real width=0.1;
real F = 1000; // Resultant of the surface force
real rho=1000;
real c1=2e6;
real c2=2e5;
real eps=5e-13;
real dt = 0.02;
real T = 1.;
int Niterations = lround(T / dt);
real NewtonTolerance=1e-7;
int NewtonIterations = 20;

// Mesh
real n=40;
real m=4;

border C0(t=0, 1){x=t*length*(n-1)/n; y=0;label=1;} // Bottom left
border C1(t=0, 1){x=length*(n-1+t)/n; y=0;label=2;} // Bottom Right
border C2(t=0, 1){x=length; y=t*width;label=2;} // Right
border C3(t=0, 1){x=(1-t/n)*length; y=width; label=2;} // Top Right
border C4(t=0, 1){x=(1-t)*(n-1)*length/n; y=width; label=1;} // Top Left
border C5(t=0, 1){x=0; y=(1-t)*width; label=3;} // Left
border C6(t=0, 1){x=(n-1)*length/n; y=width*(1-t); label=2;} // Inner Left bottom

mesh Th = buildmesh(
    C0(n-1) + C1(1) + C2(m) + C3(1) + C4(n-1) + C5(m) + C6(m)
);
//plot(C0(n-1) + C1(1) + C2(m) + C3(1) + C4(n-1) + C5(m) + C6(m), wait=true);
//plot(Th, wait=true);


// FE space
fespace Vh(Th, P1);
fespace Wh(Th, P0);

// Unknowns (on two time steps to update at each time step)
Vh dotphi1, dotphi2;
Vh olddotphi1, olddotphi2;
Vh phi1, phi2;
Vh oldphi1, oldphi2;
Wh oldp;
Wh p;

// Differentials elements for Newton's algorithm
Vh ddotphi1, ddotphi2;
Vh dphi1, dphi2;
Wh dp;


// Test functions
Vh v1, v2;
Vh w1, w2;
Wh q;

// Macros
macro Grad(v1, v2) [dx(v1), dy(v1), dx(v2), dy(v2)] //
macro cof(phi1, phi2) [dy(phi2), -dx(phi2), -dy(phi1), dx(phi1)] //

// Problem
real f = F*n/(width*length);

// Store values
real[int] NumberNewtonIterations(Niterations);
real[int] Energy(Niterations);

// First guess for the update done with Explicit Euler scheme
problem ExplicitEuler([dotphi1, dotphi2, phi1, phi2, p], [v1, v2, w1, w2, q])
    = int2d(Th)(
        rho*[dotphi1, dotphi2]' *[v1,v2]
    )
    + int2d(Th)(
        -rho*[olddotphi1, olddotphi2]' *[v1,v2]
        +dt*2*(c1+c2)*Grad(oldphi1,oldphi2)' *Grad(v1,v2)
        -dt*oldp*cof(oldphi1, oldphi2)' *Grad(v1,v2)
    )
    + int2d(Th, 0)(
        -dt*(-f*v2)
    )
    + int2d(Th)(
        [phi1,phi2]' *[w1,w2]
    )
    + int2d(Th)(
        -1*[oldphi1, oldphi2]' *[w1,w2]
        -dt*[olddotphi1,olddotphi2]' *[w1,w2]
    )
    + int2d(Th)(
        eps*p*q
    )
    + int2d(Th)(
        q*dx(oldphi1)*dy(oldphi2)
        -q*dx(oldphi2)*dy(oldphi1)
        -q
    )
    + on(3, phi1=0, phi2=0, dotphi1=0, dotphi2=0)
    ;

// Compute increment of Newton's method solving tangent problems
// i.e. has to solve F'(dotphi_{n+1}, phi_{n+1}, p_{n+1})(\delta dotphi, \delta phi, \delta p)=F(dotphi_{n+1}, phi_{n+1}, p_{n+1})
// and update (dotphi_{n+1}, phi_{n+1}, p_{n+1}) -= (\delta dotphi, \delta phi, \delta p) until ||(\delta dotphi, \delta phi, \delta p)|| < NewtonTolerance

problem ComputeNewtonUpdate([ddotphi1, ddotphi2, dphi1, dphi2, dp], [v1, v2, w1, w2, q])
    = int2d(Th)(
        rho*[ddotphi1, ddotphi2]' *[v1,v2]
        +dt*(c1+c2)*Grad(dphi1,dphi2)' *Grad(v1,v2)
        -dt*oldp*cof(dphi1,dphi2)' *Grad(v1,v2)/4
        -dt*dp*cof(oldphi1,oldphi2)' *Grad(v1,v2)/4
        -dt*dp*cof(phi1,phi2)' *Grad(v1,v2)/4
        -dt*p*cof(dphi1,dphi2)' *Grad(v1,v2)/4
    )
    + int2d(Th)(
        [dphi1, dphi2]' *[w1, w2]
        -dt*[ddotphi1, ddotphi2]' * [w1, w2]/2
    )
    + int2d(Th)(
        q*eps*dp/2
        +q*dx(oldphi1)*dy(dphi2)/4
        +q*dx(dphi1)*dy(oldphi2)/4
        +q*dx(phi1)*dy(dphi2)/4
        +q*dx(dphi1)*dy(phi2)/4
        -q*dx(oldphi2)*dy(dphi1)/4
        -q*dx(dphi2)*dy(oldphi1)/4
        -q*dx(dphi2)*dy(phi1)/4
        -q*dx(phi2)*dy(dphi1)/4
    )
    - int2d(Th)(
        rho*[dotphi1, dotphi2]' *[v1,v2]
        -rho*[olddotphi1, olddotphi2]' *[v1,v2]
        +dt*(c1+c2)*Grad(oldphi1,oldphi2)' *Grad(v1,v2)
        +dt*(c1+c2)*Grad(phi1,phi2)' *Grad(v1,v2)
        -dt*oldp*cof(oldphi1,oldphi2)' *Grad(v1,v2)/4
        -dt*oldp*cof(phi1,phi2)' *Grad(v1,v2)/4
        -dt*p*cof(oldphi1,oldphi2)' *Grad(v1,v2)/4
        -dt*p*cof(phi1,phi2)' *Grad(v1,v2)/4
    )
    - int2d(Th, 0)(
        -dt*(-f)*v2
    )
    - int2d(Th)(
        [phi1, phi2]' *[w1,w2]
        -[oldphi1, oldphi2]' *[w1,w2]
        -dt*[olddotphi1,olddotphi2]' *[w1,w2]/2
        -dt*[dotphi1,dotphi2]' *[w1,w2]/2
    )
    - int2d(Th)(
        -q
        +eps*q*oldp/2
        +eps*q*p/2
        +q*dx(oldphi1)*dy(oldphi2)/4
        +q*dx(oldphi1)*dy(phi2)/4
        +q*dx(phi1)*dy(oldphi2)/4
        +q*dx(phi1)*dy(phi2)/4
        -q*dx(oldphi2)*dy(oldphi1)/4
        -q*dx(oldphi2)*dy(phi1)/4
        -q*dx(phi2)*dy(oldphi1)/4
        -q*dx(phi2)*dy(phi1)/4
    )
    + on(3, dphi1=0, dphi2=0, ddotphi1=0, ddotphi2=0)
    ;

ofstream fenergy("output/energy.txt");
ofstream fNewtonNumber("output/NewtonNumberIterations.txt");
ofstream ferrorRelative("output/RelativeError.txt");
ofstream ferrorAbsolute("output/AbsoluteError.txt");
ofstream ferrorL2("output/L2Error.txt");
ofstream ferrorH1("output/H1Error.txt");
ofstream fdisplacement("output/vertical_displacement.txt");

ofstream fp("output/p.txt");
ofstream fphi1("output/phi1.txt");
ofstream fphi2("output/phi2.txt");
ofstream fdotphi1("output/dotphi1.txt");
ofstream fdotphi2("output/dotphi2.txt");

ofstream fdp("output/dp.txt");
ofstream fdphi1("output/dphi1.txt");
ofstream fdphi2("output/dphi2.txt");
ofstream fddotphi1("output/ddotphi1.txt");
ofstream fddotphi2("output/ddotphi2.txt");


for (int t = 0; t<Niterations; t++){

    cout << "Time: " << t*dt << "/" << T << endl << endl;

    // First guess using Explicit Euler scheme
    //ExplicitEuler;

    for (int i=0; i<NewtonIterations; i++){
        // Compute differential elements
        
        ComputeNewtonUpdate;


        // Update new configuration
        phi1[] -= dphi1[];
        phi2[] -= dphi2[];
        dotphi1[] -= ddotphi1[];
        dotphi2[] -= ddotphi2[];
        p[] -= dp[];
        
        // Compute various errors
        real errRelative = 0. ;
        real errAbsolute = 0. ;
        real errH1 = 0.;
        real errL2 = 0.;

        // Commpute relative error
        if (phi1[].linfty >0) {
            errRelative += dphi1[].linfty / phi1[].linfty;
        }
        if (phi2[].linfty >0) {
            errRelative += dphi2[].linfty / phi2[].linfty;
        }

        // Compute absolute error
        errAbsolute += dphi1[].linfty;
        errAbsolute += dphi2[].linfty;

        // Compute H1 error
        errH1 = sqrt(int2d(Th)(
            Grad(dphi1,dphi2)' *Grad(dphi1,dphi2) + dphi1*dphi1 + dphi2*dphi2
        ));

        // Compute L2 error
        errL2 = sqrt(int2d(Th)(
            dphi1*dphi1 + dphi2*dphi2
        ));
        
        // This should be more or less conserved
        real energy = int2d(Th, 0)(
                    (-f)*phi2
                )
                - int2d(Th)(
                rho*(dotphi1^2+dotphi2^2)
                +(c1+c2)*(dx(phi1)^2+dx(phi2)^2+dy(phi1)^2+dy(phi2)^2-3)
                +eps*p^2/2
                )
                ;

        // Look at the vertical displacement at the tip
        real verticalDisplacement = phi2(length, 0.5*width);

        // Write to files
        fdp << t << " " << i << " " << dp[].linfty << endl;
        fdphi1 << t << " " << i << " " << dphi1[].linfty << endl;
        fdphi2 << t << " " << i << " " << dphi2[].linfty << endl;
        fddotphi1 << t << " " << i << " " << ddotphi1[].linfty << endl;
        fddotphi2 << t << " " << i << " " << ddotphi2[].linfty << endl;
        fp << t << " " << i << " " << p[].linfty << endl;
        fphi1 << t << " " << i << " " << phi1[].linfty << endl;
        fphi2 << t << " " << i << " " << phi2[].linfty << endl;
        fdotphi1 << t << " " << i << " " << dotphi1[].linfty << endl;
        fdotphi2 << t << " " << i << " " << dotphi2[].linfty << endl;
        fenergy << t << " " << i << " " << energy << endl;
        fdisplacement << t << " " << i << " " << verticalDisplacement << endl;

        // Choose which error
        real err = errAbsolute;

        cout << "Iteration = " << i << " ; Error = " << err << " ; Tolerance = " << NewtonTolerance << " ; Vertical Displacement = " << verticalDisplacement << " ; Energy : " << energy << endl;
        if (err < NewtonTolerance){
            // Newton's iteration has converged
            // Update old time step to new time step
            cout << "Succesfully converged" << endl;

            oldphi1[] = phi1[];
            oldphi2[] = phi2[];
            olddotphi1[] = dotphi1[];
            olddotphi2[] = dotphi2[];
            oldp[] = p[];

            // If wanting to see the mesh move
            /* real coef = 1;
            real minT0 = checkmovemesh(Th, [x, y]);
            while(1){ // find a correct move mesh
                real minT = checkmovemesh(Th, [x+coef*phi1, y+coef*phi2]);
                if (minT > minT0/5) break; //if big enough
                coef /= 1.5;
            }
            mesh Thh = movemesh(Th, [x+coef*phi1, y+coef*phi2]);
            plot(Thh, wait=true); */
            plot(phi1, phi2);

            // Write how many Newton iterations it took
            fNewtonNumber << t << " " << i << endl;

            break;
        } else if (i>3 && err > 1000.){
            // Newton's iteration has diverged, early stopping
            cout << "Error exploded during Newton's iteration " << i << " at time " << t*dt << endl;
            exit(1);
        } else if (i == NewtonIterations - 1) {
            cout << "Newton's has not converged within the " << NewtonIterations << " iterations at time " << t*dt << ", error is " << err << endl;
            exit(1);
        }
    }
}

Comments on the code.

  • The mesh is constructed in two parts: the right one will be the region where the surface force is applied, while the left one is fixed on the left (Dirichlet Boundary condition for both the displacement and velocity).
  • The resultat of the force is kept constant at F=1000 N, which implies that the surface force is f=n\times F/({\rm width} \times \rm{length}), since I build the mesh so that f is applied uniformly over the width and over the last 1/n-th part of the beam, where n denotes the number of triangles along the width generated by the mesh.
  • As said in the introduction of this post, I am really, really not sure about this mesh construction. I checked that region 0 was only about the right part of the mesh, so that I only apply the force f in region 0. I hope this is not wrong.
  • For variable names, 1 respresents the x coordinate, while 2 represents the y coordinate. The prefix ‘d’ is used to denote the differential elements (the ones that should converge to 0 in Newton’s update); the prefix ‘old’ for variable names is used to store the values at iteration n, while the others (e.g. dotphi1, dotphi2, phi1, phi2, p) store the current values of the velocity/displacement/pressure during the Newton’s iterations.
  • The macro cof returns a vector representing the cofactor of the gradient of the displacement. The macro Grad also returns a vector. this is not an issue since only dot products are performed with them.
  • I wrote an ExplicitEuler scheme for a first guess for the Newton iteration, but this does not solve the problem.
  • I do not initialize any functions, since it seems they are initialized at 0 either way.
2

Can you upload the code or put on girhub? Copy from browser is a problem.
Thanks. I wanted to look at some EM problems.

3

Sorry about the delay, here it is.
midpoint.edp (9.3 KB)