0
$\begingroup$

I have quadratic finite element - shape function is quadratic. Element spans from 0 to 5.

Body force is given by (in physical coordinates) $$f_b = \int_0^5 N(x)^T b(x) dx \approx \sum_{i=1}^3 N(\xi_i)b_i(x)w_i $$ The last term is obtained from Gauss quadrature. $\xi_i$ are Gauss points and $w_i$ are corresponding weights. Say, $b=x^3$ $$ N_1(x) = ((x-x_2)(x-x_3))/( (x_1-x_2) (x_1-x_3) ) \\ N_2 (x)= ((x-x_1)(x-x_3))/( (x_2-x_1) (x_2-x_3) ) \\ N_3 (x)= ((x-x_1)(x-x_2))/( (x_3-x_1) (x_3-x_2) ) $$ When I transform the problem into isoparametric coordinates, I get

$$f_b = \sum_{i=1}^3 N(\xi_i)b_i(x)w_iJ$$

$$ N_1 = -0.5 \xi (1-\xi)\\ N_2 = (1-\xi^2)\\ N_3 = 0.5\xi(1+\xi)\\ $$

The problem is when I evaluate integrals - in physical and isoparametric coordinates - they do not match. I suspect that I am using $b(x)=x^3$ in isoparametric without converting it. But, even if I multiply $b$ by $J$, problem does not solve. I do not have much idea what is going wrong. Can someone help me?

Edit: As per comments, I checked that when $b$ is a constant, $f_b$ in both coordinates match. So problem comes when $b$ is a function of $x$. I am wondering how to convert $b$ from physical to isoparametric. To be specific, say $b=x^3$. How can I transform and integrate this in isoparametric coordinates?

| cite | improve this question | | | | |
$\endgroup$
  • 1
    $\begingroup$ Welcome to scicomp.SE! Is your code concise enough that you can post it here? My guess is that you've probably got a bug, and that it's likely in the calculation of $J$, but it's impossible to say for sure without seeing the code. Could you also put your definition of $J$ into the equations that you're using? $\endgroup$ – Tyler Olsen Jan 31 '17 at 2:38
  • $\begingroup$ I agree with Tyler, that it looks like a bug. First, try to simplify the problem. Can you compare the integrals with $b(x)=1$? If this case works fine, you have an error in values of $b$ which you use. If not, then the error is probably in the Jacobian. Then you can try to simplify the element transformation to see whether it is true. $\endgroup$ – VorKir Jan 31 '17 at 21:28
  • $\begingroup$ Edited the question. $\endgroup$ – user7423098 Feb 1 '17 at 2:18
  • $\begingroup$ In that case @Bill Greene has the answer for you. You want to find the value of $x$ that corresponds to the physical location of your quadrature point, then evaluate $b(x(\xi))$ using that value. See the answer for details. $\endgroup$ – Tyler Olsen Feb 2 '17 at 13:14
4
$\begingroup$

Instead of writing $b_i(x)$, I suggest writing it as $b(x(\xi_i))$. This is more indicative of how to evaluate your integral. $x(\xi_i)$ is the value of $x$ at the ith integration point which you can compute from the simple relationship between $x$ and $\xi$.

$$ x = \sum_{i=1}^3 N_i(\xi) x_i $$

This is the defintion of the mapping from real to parametric coordinates assumed in an isoparametric element formulation. For each integration point, $\xi_i$, you evaluate the shape functions at that point and compute the sum above to determine $x$ of the point. Then simply evaluate your $b$ function at $x$ computed from this expression.

| cite | improve this answer | | | | |
$\endgroup$
  • $\begingroup$ Can you tell me how did you arrive at the transformation $x_i = (L/2)(1+\xi_i)$ ? $\endgroup$ – user7423098 Jan 31 '17 at 23:28
  • $\begingroup$ The parametric coordinates range from -1 to +1. The real coordinates range from zero to L. So that transformation is a simple linear mapping. You want $x=0$ to match $\xi=-1$ and $x=L$ to match $\xi=1$. $\endgroup$ – Bill Greene Feb 1 '17 at 0:10
  • $\begingroup$ Okay. I had quadratic element. Also, will your transformation work for nth element? Your element starts from $x=0$ $\endgroup$ – user7423098 Feb 1 '17 at 2:14
  • $\begingroup$ In general, you have a mapping $\varphi: K \rightarrow \hat{K}$ from your element to the reference element. Then you need to take $b(x_i) = b(\varphi^{-1} (\xi_i))$ in your integral over reference element where $\xi_i$ are the integration points. This is what Bill suggests, but for particular $\varphi$, as I think. $\endgroup$ – VorKir Feb 1 '17 at 2:27
  • $\begingroup$ Sorry, no I missed the fact that you are using quadratic shape functions! I updated my answer to cover this. Note, though, that the simple expression I gave will work if the element is straight as you can easily verify by simple substitution. In the simple expression, I did assume the element started at $x=0$ so you would add $x_1$ to the result if that is not the case. $\endgroup$ – Bill Greene Feb 1 '17 at 12:20
1
$\begingroup$

An alternate approach to going from ${b}(x)$ to ${b}(x(\xi))$ by direct substitution of $x(\xi)$ in $b(x)$ is to interpolate $b(x)$ over the element using $b_i^e$ and element shape functions (quadratic) or even linear shape functions. So, one would write $$ b(x(\xi)) \approx \sum_{i=1}^3 b_i^e N_i(\xi) \text{ on } \Omega^e $$ with $b_i^e$ being the nodal value of the body force at node i on element $e$. Of course, when $b(x)$ is constant or linear, this would be exact.

| cite | improve this answer | | | | |
$\endgroup$
  • $\begingroup$ It seems like the only thing that could happen is that you lose accuracy by doing this. You'll never gain any information, even if your shape functions are of a higher polynomial degree than the body force. In the case that your shape functions are of lower order, you're approximating the body force with a lower-order piecewise polynomial. This method would probably converge as you refine the mesh, but it seems like there's no point in deliberately reducing the accuracy if you can avoid it, considering that it's practically zero effort to do the $b(x(\xi))$ method. $\endgroup$ – Tyler Olsen Feb 5 '17 at 13:57
  • $\begingroup$ The integration being carried out is numerical (Gauss quadrature). With Gauss quadrature, an $n$-point quadrature integrates a polynomial of order $2n-1$ at most, exactly. The OP is using the three-point quadrature rule which means a polynomial of at most a fifth order can be integrated exactly in general. Now, let's take a look at the integrand: $N(\xi)*x^\prime(\xi)$ is a polynomial of cubic order. Therefore, the function $b(\xi)$ can at most be quadratic for exact integration. $\endgroup$ – hpc Feb 5 '17 at 15:55

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.