On Tue, Jan 23, 2018 at 11:14 AM, Yann Jobic <yann.jo...@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/>