3
$\begingroup$

I'm trying to develop a kind of nonlinear FEM application using mathematica to solve a bvp like the following:

$$ \gamma(u') ~u^{iv} + 2 \gamma'(u') u''' u''+ u''^3 = f(x) $$ where $u = \tilde{u}(x)$, $\gamma = \tilde{\gamma}(u')$ and $'$ denotes derivative with respect to the function variable, thus:

$$ u' = \frac{\text{d} u(x)}{\text{d}x} \quad \quad \gamma' = \frac{\text{d} \gamma(u')}{\text{d} u'} $$

The next step is getting the weak form:

$$ \int_{\Omega} \left( \gamma ~u^{iv} + 2 \gamma' u''' u''+ u''^3 \right) \psi = \int_{\Omega} f(x)\psi $$

where $\Omega$ is the domain and $\psi = \tilde{\psi}(x)$ is the test function. In my code, I'm using a Cubic Hermite Shape Functions, then the approximate function reads as follows:

$$ \hat{u} = \sum_i \alpha_i \psi_i^{(a)} + \beta_i \psi_i^{(b)} $$

where $\psi_i^{(a)}$ and $\psi_i^{(b)}$ are piecewise functions, and $\alpha_i$ and $\beta_i$ are coefficients related to the values of function and its first derivative, respectively. Let bear in mind the chain rule:

$$ \frac{\text{d} \gamma(u')}{\text{d} x} = \frac{\text{d} \gamma(u')}{\text{d} u'} \frac{\text{d} u'}{\text{d} x} = \gamma' u'' $$

After some calculation, the weak form assumes the following form:

$$ \int_{\Omega}\gamma~u^{iv} = \left[\left(\gamma u''' + \frac{1}{2}\gamma' u''^2 \right)\psi-\left(\gamma \psi \right)' u'' \right]_{\partial \Omega} + \int_{\Omega} \frac{1}{2} \gamma'' u''^3 \psi + \frac{3}{2} \gamma' u''^2 \psi' + \gamma u'' \psi '' $$

$$ \int_{\Omega} 2 \gamma' u''' u'' \psi = \left[\gamma' u''^2 \psi \right]_{\partial \Omega} - \int_{\Omega} \gamma'' \psi u''^3 \psi + \gamma' u''^2 \psi' $$

As starting point, however, I'm using a linear equation, setting $\gamma(u') = 1 $ (but at the end, the procedure is the same). Then, the simple equation is:

$$ u^{iv} = f(x) \Longrightarrow \int_{\Omega} u^{iv} \psi = \int_{\Omega} f(x) \psi $$ then $$ \left[u''' \psi - u'' \psi' \right]_{\partial \Omega} + \int_{\Omega} u'' \psi'' - \int_{\Omega} f(x) \psi = 0 $$ If the domain is discretised in n points (n-1 elements), this equation represents a system of 2n (non)linear equations. Let suppose to apply BCs on the function and its first derivative at both edges. The system reads as follows:

$$ R(1) = \alpha_1 - u(x_0) \\\\ R(2) = \beta_1 - u'(x_0) \\\\ ... \\\\ R(i) = \int_{\Omega_i} u'' \psi'' - \int_{\Omega_i} f(x) \psi \\\\ ... \\\\ R(2n-1) = \alpha_n - u(x_1) \\\\ R(2n) = \beta_n - u'(x_1) \\\\ $$

where the domain is $\Omega = [x_0, x_1]$. Now the goal is solving this system (in general nonlinear). I'm trying to do that using FindRoot. Thus I defined the following functions:

(* integration of forcing term *)
funFF[X_, f_] := Module[
  {F, n},
  n = Length[X];
  F = ConstantArray[0, 2 n];
  Do[
    jr = 2 j - 1;
    F[[jr ;; jr + 3]] = 
    F[[jr ;; jr + 3]] + 
    Flatten[NIntegrate[
        f[x] NN0[x, X[[j]], X[[j + 1]]], {x, X[[j]], X[[j + 1]]}]],
    {j, 1, n - 1}
    ];
  F
  ];

  (* assembly of system *)
  SYS1[X_, U_?VectorQ, RHS_, fun_] := Module[
    {EQ, x1, x2, f1, p1, f2, p2, F, P, R, n},
    n = Length[X]; 
    (* value of function *)
    F = U[[1 ;; Length[U] ;; 2]];
    (* value of first derivative *)
    P = U[[2 ;; Length[U] ;; 2]];
    (* initialitazion *)
    EQ = ConstantArray[0, 2 n];
    (* starting loop *)
    Do[
      jr = 2 j - 1;
      (* values for the current element *)
      x1 = X[[j]];
      x2 = X[[j + 1]];
      If[Length[F] == 0,
        f1 = F;
        f2 = F,
        (* if non-constant *)
        f1 = F[[j]];
        f2 = F[[j + 1]]
      ];
      If[Length[P] == 0,
        p1 = P;
        p2 = P,
        (* if non-constant *)
        p1 = P[[j]];
        p2 = P[[j + 1]]
      ];
      (* j-th equation *)
      EQ[[jr ;; jr + 3]] = EQ[[jr ;; jr + 3]] + 
             Flatten[NIntegrate[fun[x, x1, x2, f1, p1, f2, p2], {x, x1, x2}]],
      {j, 1, n - 1}
    ];
    (* assembly of system *)
    R = EQ - RHS;
    (* boundary conditions *)
    R[[1]] = U[[1]] - wf0;
    R[[2]] = U[[2]] - wp0;
    R[[2 n - 1]] = U[[2 n - 1]] - wf1;
    R[[2 n]] = U[[2 n]] - wp1;
    (*output*)
    R
];

where I assumed that $w(x)$ is the exact solution, thus wf0 = w(x0), wp0 = w'(x0), wf1 = w(x1)$ and wp1 = w'(x1). Then:

(* forcing term *)
f1[x_] = w''''[x];

(* integration of forcing term *)
F1 = funFF[XX, f1];

(* function to be integrated*)
fun1[x_, x1_, x2_, f1_, p1_, f2_, p2_] := Ne2[x, x1, x2, f1, p1, f2, p2] NN2[x, x1, x2];

(* solution *)
U1 = UU /. Flatten[FindRoot[SYS1[XX, UU, F1, fun1] == 0, {UU, U0}]

where:

  • XX is the grid
  • U0 = {\alpha_1, \beta_1, \alpha_2, \beta_2, ... } is the guess
  • Ne2 contains the interpolation inside an element
  • NN2 contains the four shape functions
  • fun is the integrand for each element, i.e $u'' \psi ''$

This simple example works, but it is really slow even if it's a linear case. How can I improve/optimize the code? I'm new in this kind of programming, so I'm sure I'm missing some speeding-up tricks.

Any suggestion is very appreciated.

Best,

Petrus

$\endgroup$
  • 1
    $\begingroup$ In the finite element method the numerical integration over the element is typically done using very minimal sampling (i.e. one or two quadrature points). You should not be using NIntegrate , rather work out by hand the 1- or 2- point gauss integration formula and use that. $\endgroup$ – george2079 Jul 11 '14 at 20:39
  • $\begingroup$ Having taken @george2079's suggestion, your code should be ideal for Compile, which will speed it up significantly. $\endgroup$ – Oleksandr R. Dec 26 '14 at 14:16
  • $\begingroup$ @Petrus, where does your project stand? Do you now have a reliable, reasonably fast procedure? $\endgroup$ – bbgodfrey Mar 1 '15 at 15:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.