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
