<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
  </head>
  <body text="#000000" bgcolor="#FFFFFF">
    <p>Hello Matt,<br>
    </p>
    <br>
    <div class="moz-cite-prefix">Le 15/02/2018 à 00:41, Matthew Knepley
      a écrit :<br>
    </div>
    <blockquote type="cite"
cite="mid:CAMYG4GnD_7mbx5WHoJ5raQrSr_sj3VZs-EvsaP4f_jWV-OkGEg@mail.gmail.com">
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
      <div dir="ltr">
        <div class="gmail_extra">
          <div class="gmail_quote">On Tue, Jan 23, 2018 at 11:14 AM,
            Yann Jobic <span dir="ltr"><<a
                href="mailto:yann.jobic@univ-amu.fr" target="_blank"
                moz-do-not-send="true">yann.jobic@univ-amu.fr</a>></span>
            wrote:<br>
            <blockquote class="gmail_quote" style="margin:0 0 0
              .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello,<br>
              <br>
              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 </blockquote>
            <div><br>
            </div>
            <div>I need to understand what you mean by "quadrature
              points". I mean the following thing:</div>
            <div><br>
            </div>
            <div>  I want to do an integral over the domain for a
              variational form:</div>
            <div>    <v, f(u)> = \int_\Omega v . f(u, x)</div>
            <div>  Now I can break this up into a sum of integrals over
              each element because integrals are additive</div>
            <div>    = \sum_T \int_T v . f(u)</div>
            <div>  And we normally integrate over a reference element
              T_r instead</div>
            <div>    = \sum_T \int_{T_r} v_r . f_r(u_r, x) |J|</div>
            <div>  And then we approximate these cell integrals with
              quadrature</div>
            <div>    = \sum_T \sum_q v_r(x_q) . f_r(u_r(x_q), x_q)
              |J(x_q)| w_q</div>
            <div>  The quadrature points x_q and weights w_q are defined
              on the reference element. This means they</div>
            <div>  are not shared by definition.</div>
            <div><br>
            </div>
            <div>Does this make sense?</div>
          </div>
        </div>
      </div>
    </blockquote>
    Perfectly. I was mixing the notions of nodes of an element, and the
    quadrature points. Thanks a lot for the clarification !<br>
    <br>
    Best regards,<br>
    <br>
    Yann<br>
    <blockquote type="cite"
cite="mid:CAMYG4GnD_7mbx5WHoJ5raQrSr_sj3VZs-EvsaP4f_jWV-OkGEg@mail.gmail.com">
      <div dir="ltr">
        <div class="gmail_extra">
          <div class="gmail_quote">
            <div><br>
            </div>
            <div>  Thanks,</div>
            <div><br>
            </div>
            <div>     Matt</div>
            <div> </div>
            <blockquote class="gmail_quote" style="margin:0 0 0
              .8ex;border-left:1px #ccc solid;padding-left:1ex">i
              understood so far (which is not far...)<br>
              I took the example of the jacobian calculus.<br>
              <br>
              I found this comment in dmplexsnes.c, which explains the
              basic idea:<br>
              1725:   /* 1: Get sizes from dm and dmAux */<br>
              1726:   /* 2: Get geometric data */<br>
              1727:   /* 3: Handle boundary values */<br>
              1728:   /* 4: Loop over domain */<br>
              1729:   /*   Extract coefficients */<br>
              1730:   /* Loop over fields */<br>
              1731:   /*   Set tiling for FE*/<br>
              1732:   /*   Integrate FE residual to get elemVec */<br>
              [...]<br>
              1740:   /* Loop over domain */<br>
              1741:   /*   Add elemVec to locX */<br>
              <br>
              I almost get that. The critical part should be :<br>
              loop over fieldI<br>
              2434:     PetscFEGetQuadrature(fe, &quad);<br>
              2435:     PetscFEGetDimension(fe, &Nb);<br>
              2436:     PetscFEGetTileSizes(fe, NULL, &numBlocks,
              NULL, &numBatches);<br>
              2437:     PetscQuadratureGetData(quad, NULL, NULL,
              &numQuadPoints, NULL, NULL);<br>
              2438:     blockSize = Nb*numQuadPoints;<br>
              2439:     batchSize = numBlocks * blockSize;<br>
              2440:     PetscFESetTileSizes(fe, blockSize, numBlocks,
              batchSize, numBatches);<br>
              2441:     numChunks = numCells / (numBatches*batchSize);<br>
              2442:     Ne        = numChunks*numBatches*batchSize<wbr>;<br>
              2443:     Nr        = numCells % (numBatches*batchSize);<br>
              2444:     offset    = numCells - Nr;<br>
              2445:     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {<br>
              <br>
              From there, we can have the numbering with (in dtfe.c)<br>
              basic idea :<br>
              6560: $   Loop over element matrix entries (f,fc,g,gc
              --> i,j):<br>
              Which leads to :<br>
              4511:       PetscPrintf(PETSC_COMM_SELF, "Element matrix
              for fields %d and %d\n", fieldI, fieldJ);<br>
              4512:       for (fc = 0; fc < NcI; ++fc) {<br>
              4513:         for (f = 0; f < NbI; ++f) {<br>
              4514:           const PetscInt i = offsetI + f*NcI+fc;<br>
              4515:           for (gc = 0; gc < NcJ; ++gc) {<br>
              4516:             for (g = 0; g < NbJ; ++g) {<br>
              4517:               const PetscInt j = offsetJ + g*NcJ+gc;<br>
              4518:               PetscPrintf(PETSC_COMM_SELF, "
              elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc,
              PetscRealPart(elemMat[eOffset+<wbr>i*totDim+j]));<br>
              [...]<br>
              4525:     cOffset    += totDim;<br>
              4526:     cOffsetAux += totDimAux;<br>
              4527:     eOffset    += PetscSqr(totDim);<br>
              4528:   }<br>
              <br>
              But i didn't get how you can find that there are
              duplicates quadrature nodes, and how you manage them.<br>
              Maybe i looked at the wrong part of the code ?<br>
              <br>
              Thanks !<br>
              <br>
              Best regards,<br>
              <br>
              Yann<br>
              <br>
              <br>
              ---<br>
              L'absence de virus dans ce courrier électronique a été
              vérifiée par le logiciel antivirus Avast.<br>
              <a href="https://www.avast.com/antivirus" rel="noreferrer"
                target="_blank" moz-do-not-send="true">https://www.avast.com/antiviru<wbr>s</a><br>
              <br>
            </blockquote>
          </div>
          <br>
          <br clear="all">
          <div><br>
          </div>
          -- <br>
          <div class="gmail_signature" data-smartmail="gmail_signature">
            <div dir="ltr">
              <div>
                <div dir="ltr">
                  <div>What most experimenters take for granted before
                    they begin their experiments is infinitely more
                    interesting than any results to which their
                    experiments lead.<br>
                    -- Norbert Wiener</div>
                  <div><br>
                  </div>
                  <div><a href="http://www.caam.rice.edu/%7Emk51/"
                      target="_blank" moz-do-not-send="true">https://www.cse.buffalo.edu/~knepley/</a><br>
                  </div>
                </div>
              </div>
            </div>
          </div>
        </div>
      </div>
    </blockquote>
    <br>
    <pre class="moz-signature" cols="72">-- 
___________________________

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
</pre>
  <div id="DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2"><br /> <table style="border-top: 1px solid #D3D4DE;">
        <tr>
      <td style="width: 55px; padding-top: 18px;"><a href="https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient" target="_blank"><img src="https://ipmcdn.avast.com/images/icons/icon-envelope-tick-round-orange-animated-no-repeat-v1.gif" alt="" width="46" height="29" style="width: 46px; height: 29px;" /></a></td>
                <td style="width: 470px; padding-top: 17px; color: #41424e; font-size: 13px; font-family: Arial, Helvetica, sans-serif; line-height: 18px;">Garanti sans virus. <a href="https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient" target="_blank" style="color: #4453ea;">www.avast.com</a>               </td>
        </tr>
</table>
<a href="#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2" width="1" height="1"> </a></div></body>
</html>