[petsc-users] FEM & conformal mesh

[Matthew Knepley]

On Tue, Jan 23, 2018 at 11:14 AM, Yann Jobic <yann.jobic at univ-amu.fr> wrote:

> Hello,
>
> 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


I need to understand what you mean by "quadrature points". I mean the
following thing:

  I want to do an integral over the domain for a variational form:
    <v, f(u)> = \int_\Omega v . f(u, x)
  Now I can break this up into a sum of integrals over each element because
integrals are additive
    = \sum_T \int_T v . f(u)
  And we normally integrate over a reference element T_r instead
    = \sum_T \int_{T_r} v_r . f_r(u_r, x) |J|
  And then we approximate these cell integrals with quadrature
    = \sum_T \sum_q v_r(x_q) . f_r(u_r(x_q), x_q) |J(x_q)| w_q
  The quadrature points x_q and weights w_q are defined on the reference
element. This means they
  are not shared by definition.

Does this make sense?

  Thanks,

     Matt


> i understood so far (which is not far...)
> I took the example of the jacobian calculus.
>
> I found this comment in dmplexsnes.c, which explains the basic idea:
> 1725:   /* 1: Get sizes from dm and dmAux */
> 1726:   /* 2: Get geometric data */
> 1727:   /* 3: Handle boundary values */
> 1728:   /* 4: Loop over domain */
> 1729:   /*   Extract coefficients */
> 1730:   /* Loop over fields */
> 1731:   /*   Set tiling for FE*/
> 1732:   /*   Integrate FE residual to get elemVec */
> [...]
> 1740:   /* Loop over domain */
> 1741:   /*   Add elemVec to locX */
>
> I almost get that. The critical part should be :
> loop over fieldI
> 2434:     PetscFEGetQuadrature(fe, &quad);
> 2435:     PetscFEGetDimension(fe, &Nb);
> 2436:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
> 2437:     PetscQuadratureGetData(quad, NULL, NULL, &numQuadPoints, NULL,
> NULL);
> 2438:     blockSize = Nb*numQuadPoints;
> 2439:     batchSize = numBlocks * blockSize;
> 2440:     PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize,
> numBatches);
> 2441:     numChunks = numCells / (numBatches*batchSize);
> 2442:     Ne        = numChunks*numBatches*batchSize;
> 2443:     Nr        = numCells % (numBatches*batchSize);
> 2444:     offset    = numCells - Nr;
> 2445:     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
>
> From there, we can have the numbering with (in dtfe.c)
> basic idea :
> 6560: $   Loop over element matrix entries (f,fc,g,gc --> i,j):
> Which leads to :
> 4511:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and
> %d\n", fieldI, fieldJ);
> 4512:       for (fc = 0; fc < NcI; ++fc) {
> 4513:         for (f = 0; f < NbI; ++f) {
> 4514:           const PetscInt i = offsetI + f*NcI+fc;
> 4515:           for (gc = 0; gc < NcJ; ++gc) {
> 4516:             for (g = 0; g < NbJ; ++g) {
> 4517:               const PetscInt j = offsetJ + g*NcJ+gc;
> 4518:               PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]:
> %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
> [...]
> 4525:     cOffset    += totDim;
> 4526:     cOffsetAux += totDimAux;
> 4527:     eOffset    += PetscSqr(totDim);
> 4528:   }
>
> But i didn't get how you can find that there are duplicates quadrature
> nodes, and how you manage them.
> Maybe i looked at the wrong part of the code ?
>
> Thanks !
>
> Best regards,
>
> Yann
>
>
> ---
> L'absence de virus dans ce courrier électronique a été vérifiée par le
> logiciel antivirus Avast.
> https://www.avast.com/antivirus
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180214/7135064d/attachment.html>