<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Tue, Jan 23, 2018 at 11:14 AM, Yann Jobic <span dir="ltr"><<a href="mailto:yann.jobic@univ-amu.fr" target="_blank">yann.jobic@univ-amu.fr</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello,<br>
<br>
I'm trying to understand the numbering of quadrature points in order to solve the FEM system, and how you manage this numbering in order to allow conformal mesh. I looked in several files in order to understand. Here's what </blockquote><div><br></div><div>I need to understand what you mean by "quadrature points". I mean the following thing:</div><div><br></div><div>  I want to do an integral over the domain for a variational form:</div><div>    <v, f(u)> = \int_\Omega v . f(u, x)</div><div>  Now I can break this up into a sum of integrals over each element because integrals are additive</div><div>    = \sum_T \int_T v . f(u)</div><div>  And we normally integrate over a reference element T_r instead</div><div>    = \sum_T \int_{T_r} v_r . f_r(u_r, x) |J|</div><div>  And then we approximate these cell integrals with quadrature</div><div>    = \sum_T \sum_q v_r(x_q) . f_r(u_r(x_q), x_q) |J(x_q)| w_q</div><div>  The quadrature points x_q and weights w_q are defined on the reference element. This means they</div><div>  are not shared by definition.</div><div><br></div><div>Does this make sense?</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">i understood so far (which is not far...)<br>
I took the example of the jacobian calculus.<br>
<br>
I found this comment in dmplexsnes.c, which explains the basic idea:<br>
1725:   /* 1: Get sizes from dm and dmAux */<br>
1726:   /* 2: Get geometric data */<br>
1727:   /* 3: Handle boundary values */<br>
1728:   /* 4: Loop over domain */<br>
1729:   /*   Extract coefficients */<br>
1730:   /* Loop over fields */<br>
1731:   /*   Set tiling for FE*/<br>
1732:   /*   Integrate FE residual to get elemVec */<br>
[...]<br>
1740:   /* Loop over domain */<br>
1741:   /*   Add elemVec to locX */<br>
<br>
I almost get that. The critical part should be :<br>
loop over fieldI<br>
2434:     PetscFEGetQuadrature(fe, &quad);<br>
2435:     PetscFEGetDimension(fe, &Nb);<br>
2436:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);<br>
2437:     PetscQuadratureGetData(quad, NULL, NULL, &numQuadPoints, NULL, NULL);<br>
2438:     blockSize = Nb*numQuadPoints;<br>
2439:     batchSize = numBlocks * blockSize;<br>
2440:     PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);<br>
2441:     numChunks = numCells / (numBatches*batchSize);<br>
2442:     Ne        = numChunks*numBatches*batchSize<wbr>;<br>
2443:     Nr        = numCells % (numBatches*batchSize);<br>
2444:     offset    = numCells - Nr;<br>
2445:     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {<br>
<br>
>From there, we can have the numbering with (in dtfe.c)<br>
basic idea :<br>
6560: $   Loop over element matrix entries (f,fc,g,gc --> i,j):<br>
Which leads to :<br>
4511:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);<br>
4512:       for (fc = 0; fc < NcI; ++fc) {<br>
4513:         for (f = 0; f < NbI; ++f) {<br>
4514:           const PetscInt i = offsetI + f*NcI+fc;<br>
4515:           for (gc = 0; gc < NcJ; ++gc) {<br>
4516:             for (g = 0; g < NbJ; ++g) {<br>
4517:               const PetscInt j = offsetJ + g*NcJ+gc;<br>
4518:               PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+<wbr>i*totDim+j]));<br>
[...]<br>
4525:     cOffset    += totDim;<br>
4526:     cOffsetAux += totDimAux;<br>
4527:     eOffset    += PetscSqr(totDim);<br>
4528:   }<br>
<br>
But i didn't get how you can find that there are duplicates quadrature nodes, and how you manage them.<br>
Maybe i looked at the wrong part of the code ?<br>
<br>
Thanks !<br>
<br>
Best regards,<br>
<br>
Yann<br>
<br>
<br>
---<br>
L'absence de virus dans ce courrier électronique a été vérifiée par le logiciel antivirus Avast.<br>
<a href="https://www.avast.com/antivirus" rel="noreferrer" target="_blank">https://www.avast.com/antiviru<wbr>s</a><br>
<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.caam.rice.edu/~mk51/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div>
</div></div>