[petsc-users] FEM & conformal mesh

[Yann Jobic]

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