Hello Matt,
Le 15/02/2018 à 00:41, Matthew Knepley a écrit :
On Tue, Jan 23, 2018 at 11:14 AM, Yann Jobic <[email protected]
<mailto:[email protected]>> 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
