Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
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. –  george2079 Jul 11 at 20:39

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.