3
$\begingroup$

I am implementing a MatLab program to solve the equation given in this paper, which involves solving integrals coming from the variational formulation of the problem. One of them is $$\int_{\partial \Omega} f\left(\sum_{i=1}^Nu_i\phi_i\right) \cdot \phi_j d\Omega$$ where $\phi_j$ is a function such that, in the node $(x_j, y_j)$ of the triangular partition of $\Omega$, its value is 1 and 0 elsewhere. That is, $\{\phi_i\}$ is the standard nodal basis defined for a triangular partition on FEM.

My confusion is how to evaluate this integral. Here is what I think I need to do:

At each iteration, we want to evaluate $u^{(k)} = (u_1^{(k)}, ..., u_N^{(k)})^T$ an approximation of $u$ at each node from our triangular partition with $N$ nodes. This integral is only evaluated at the nodes which lie on the boundary $\partial \Omega$. Suppose, for example, that $(x_1, y_1)$ is one of those nodes and we want the value of $u$ at this node (ie, we want u_1 at some iteration). Then, we should have to evaluate:

$$\int_{\partial \Omega} f\left(\sum_{i=1}^Nu_i\phi_i(x_1, y_1)\right) \cdot \phi_1(x_1, y_1) d\Omega = \int_{\partial \Omega} f(u_1) d\Omega$$

I am really confused if this is correct, and if it is, how do I calculate this integral?

I am sorry if the question is confusing, but the paper can give a clear idea of what should be done, if needed.

Thanks!

$\endgroup$

2 Answers 2

3
$\begingroup$

So taking your first expression, and giving it a name to make an equation out of it, and substituting a quadrature rule says: $$ \big[F(u)\big]_{j} \approx \sum_q w_q f\Big(\sum_i u_i \phi_i(x_q)\Big)\phi_j(x_q) $$ where $w_q$ and $x_q$ are the weights and points of your quadrature rule. If this is the RHS of your overall system, then $\phi_j$ should probably not appear since this would make a matrix for all $i$ and $j$ in your domain. If this is in the typical left-hand side, then you have additional contribution to the $i,j$ entry of your stiffness matrix. If $x_q$ is not in any of the elements surrounding $x_i$, then you can skip that term in the inner sum since $\phi_i(x_q)=0$ in that case and just add $f(0)$ to your partial sum for that $i$. It might be good if you edited to include the whole equation you're trying to solve, e.g. $Lu=F(u)$, for a linear operator $L$, or $F(u)={\rm something else}$ to help clarify. Either way, you're probably looking to form a matrix equation to solve for each $u_i$. If $u_i$ is unknown and $f$ is non-linear, then forming a matrix out of this system to solve for $u_i$ is probably impossible.

Due to this inner sum inside $f$ which is presumably non-linear, this is going to be an expensive term to compute. You might be able to quickly build the full $N_qN_i$ table if both $N$s are of reasonable size and $f$ is cheap.

$\endgroup$
2
$\begingroup$

You seem to be confused with the definition of the basis $\{\phi_i\}.$ You are right that $\phi_i(x_i,y_i)=1$ and $\phi_i(x_j,y_j)=0$ for all $i\neq j.$ What you probably miss is that $\phi_i$ is linear, i.e. it is non-zero on all triangular elements which touch $(x_i,y_i).$

The paper describes a standard Galerkin discretisation method. The integral in the left-hand side of your expression can be evaluated analytically, since the function $f(u)=-u(u-a)(u-b)$ is a polynomial.

$\endgroup$
6
  • $\begingroup$ Thanks for your answer. Indeed, I forgot about this property of nodal basis, but I am still confused about how to evaluate $F(u)$, using that $u$ is approximation of the sum of $u_i's$. Can you indicate a book or paper that gives a brief explain of the method or an example of how to calculate a integral like that? $\endgroup$
    – Giiovanna
    Nov 2, 2015 at 8:00
  • $\begingroup$ I think my confusion is not clear: In FEM, we usually define a local basis, ie, we define a basis for each triangle using a determinant formula, for example. Suppose I generate this base. How can I evaluate this integral, since it seems to use a notion of "global basis"? $\endgroup$
    – Giiovanna
    Nov 2, 2015 at 11:26
  • $\begingroup$ Each basis function is a function of position, so you can evaluate it anywhere you like within your domain including on the boundary. What's typically done is to set up a collection of quadrature points on the boundary, on each face or edge on the boundary as appropriate. Then you simply integrate as you would in the full domain. $\endgroup$
    – Bill Barth
    Nov 2, 2015 at 15:46
  • $\begingroup$ @BillBarth I am using the idea of local basis function, as given here:books.google.com.br/… . My $u_i$ must be the approximate value of $u$ at node $(x_i, y_i)$, right? So, aren't the nodes the points that I should use? $\endgroup$
    – Giiovanna
    Nov 2, 2015 at 20:58
  • $\begingroup$ @Giiovanna, I can't read that link. Google forbids me, and I don't read Portuguese. Yes, but unless you are using a nodal integration rule, you will also have integration points in the interior of each element or the interior of each face that you can conduct integration over in the usual FEM way though one dimension lower. If you replace your continuous integral with a quadrature rule, you should get some familiar looking matrix entries for boundary nodes. $\endgroup$
    – Bill Barth
    Nov 2, 2015 at 21:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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