[petsc-users] FEM & conformal mesh

[Yann Jobic]

Hello Matt,


Le 15/02/2018 à 00:41, Matthew Knepley a écrit :
> On Tue, Jan 23, 2018 at 11:14 AM, Yann Jobic <yann.jobic at univ-amu.fr 
> <mailto: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?
Perfectly. I was mixing the notions of nodes of an element, and the 
quadrature points. Thanks a lot for the clarification !

Best regards,

Yann
>
>   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 <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/%7Emk51/>

-- 
___________________________

Yann JOBIC
HPC engineer
IUSTI-CNRS UMR 7343 - Polytech Marseille
Technopôle de Château Gombert
5 rue Enrico Fermi
13453 Marseille cedex 13
Tel : (33) 4 91 10 69 43
Fax : (33) 4 91 10 69 69



---
L'absence de virus dans ce courrier électronique a été vérifiée par le logiciel antivirus Avast.
https://www.avast.com/antivirus
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180215/cd5b79b4/attachment.html>