Higher-order elements: how is this supposed to work?

mennodeij commented 1 year ago

I'm trying to use the higher-order elements functionality but I'm getting all kind of strange errors (floating point errors, index out-of-range, etc.). There is no example for how to use this, can you maybe clarify which calls should be made on the tiogaInterface for this to work? (or if this functionality is not supported anymore, knowing that'd help me move on)

jsitaraman commented 1 year ago

To my knowledge, this should still work. Let me look for an example and get back to you.

On Tue, Feb 7, 2023 at 5:20 AM Menno Deij - van Rijswijk < @.***> wrote:

I'm trying to use the higher-order elements functionality but I'm getting all kind of strange errors (floating point errors, index out-of-range, etc.). There is no example for how to use this, can you maybe clarify which calls should be made on the tiogaInterface for this to work? (or if this functionality is not supported anymore, knowing that'd help me move on)

— Reply to this email directly, view it on GitHub https://github.com/jsitaraman/tioga/issues/67, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACH3YMDJXBJHNW6MUQLRXKDWWJDY3ANCNFSM6AAAAAAUT6U6C4 . You are receiving this because you are subscribed to this thread.Message ID: @.***>

mennodeij commented 1 year ago

Ok, thank you that would be very helpful 🙂

mennodeij commented 1 year ago

I think I understand what is going on. A higher-order element can be defined, as long as it has at most 8 vertices (hexa). I was not working on that assumption.

ackirby88 commented 1 year ago

@mennodeij For high order elements, you can have (almost) arbitrarily higher-order elements for the solution array (p=10, i.e., up to 11 points in each coordinate direction).

In order to utilize the high order functionality, you need to do roughly the following steps:

[ ] Register the grid data with TIOGA by calling tioga_registergrid_data_:

void tioga_registergrid_data_(int *btag,int *nnodes,double *xyz,int *iblank,int *nwbc, int *nobc,int *wbcnode,
              int *obcnode,int *ntypes,...);

btag: body tag index (needs to be globally unique per body) nnodes: total number of unique simplex mesh nodes on this mpi process (note: do not include high-order curved geometry nodes -- that comes later) xyz: physical coordinate array for the nodes, stacked as (x1,y1,z1,x2,y2,z2,...), size=3*nnodes iblank: iblanking array for the nodes, size=nnodes nwbc: number of wall boundary condition nodes nobc: number of overset boundary condition nodes wbcnode: wall boundary condition node indices, size=nwbc obcnode: overset boundary condition node indices, size=nobc ntypes: number of element types used (see example below) ...: lists of element stride, number of each element, element node connectivity

Example for both low-order and high-order codes:

  /* every proc has to register an unstructured grid even if there is nothing in it */
  d->tioga->registergrid_data(&d->global_body_tag,
                              &f->nnode,
                               f->xgeom_translated,
                               f->iblank,
                              &f->nwbc,
                              &f->nobc,
                               f->iwbcnode,
                               f->iobcnode,
                              &f->ncell_types, // 4 (tets/pyramids/prisms/hexes)
                              &f->kstride4,    // tet striding = 4
                              &f->ntet,        // number of rank-local tets
                               f->ndc4,        // tet connectivity (node indices composing each tet)
                              &f->kstride5,    // pyramid striding = 5
                              &f->npyramid,    // number of rank-local pyramids
                               f->ndc5,        // pyramid connectivity (node indices composing each pyramid)
                              &f->kstride6,    // prism striding = 6
                              &f->nprism,      // number of rank-local prisms
                               f->ndc6,        // prism connectivity (node indices composing each prism)
                              &f->kstride8,    // hex striding = 8
                              &f->nhex,        // number of rank-local hexes
                               f->ndc8);       // hex connectivity (node indices composing each hex)

[ ] Tell TIOGA that you have a high-order iblank_cell array by calling the function void tioga_setcelliblank_

void tioga_setcelliblank_(int *iblank_cell); // iblank_cell has length equal total number of mesh elements on mpi rank

[ ] Tell TIOGA that you have high-order data by calling the function tioga_set_highorder_callback_:

void tioga_set_highorder_callback_(count_receptor_nodes,
                                   create_receptor_nodes,
                                   donor_inclusion_test,
                                   create_donor_frac,
                                   convert_to_receptor_coefficients);

where the function inputs are five callback functions:

void count_receptor_nodes(int *cell_index,int *num_receptor_nodes);

/** Driver interface callback function:
*      callback function for Tioga that counts receptor nodes
*
* @param [in]  cell_index          cell index local to processor, note base 1
* @param [out] num_receptor_nodes  number of receptor nodes for cell with index \a cell_index, (e.g., p=3 hex is 64 nodes)
*/
void count_receptor_nodes(int *cell_index,int *num_receptor_nodes);

void create_receptor_nodes(int *cell_index,int *num_receptor_nodes,double *receptor_nodes);

/** Driver interface callback function:
*      callback function for Tioga that creates receptor nodes
*
* @param [in]  cell_index          cell index local to processor, base 1
* @param [in]  num_receptor_nodes  number of receptor nodes for cell with index \a cell_index can be used to allocate
* @param [out] receptor_nodes      location of receptor nodes 3*num_receptor_nodes: e.g., [x1,y1,z1,x2,y2,z2,...,xN,yN,zN] where N=num_receptor_nodes
*/
void create_receptor_nodes(int *cell_index,int *num_receptor_nodes,double *receptor_nodes);

void donor_inclusion_test(int *cell_index,double *pt_xyz,int *pass_flag,Real *rst_out);

/** Driver interface callback function:
*      callback function for Tioga that tests high order inclusion of a point with a cell
*
* @param [in]  cell_index  index of the cell testing for inclusion
* @param [in]  pt_xyz      physical location of point used for inclusion (length 3)
* @param [out] pass_flag   set to 1 if \a xyz is inside cell and 0 if it is not
* @param [out] rst_out     if \a pass_flag is 1 then put the natural coordinates here (-1<r,s,t<1, length 3)
*/
void donor_inclusion_test(int *cell_index,double *pt_xyz,int *pass_flag,Real *rst_out);

void create_donor_frac(int *cell_index,double *xyz,int *nfrac,int *index,double *frac,double *rst_in,int *dimfrac);

/** Driver interface callback function:
*      callback function for Tioga that creates interpolation weights
*
* @param [in]  cell_index  index of cell creating interpolation weights for
* @param [in]  xyz         physical location (not used anymore since \a rst_in is always passed in)
* @param [out] nfrac       number of modes in the basis used for interpolation
* @param [out] index       index to access solution array index is base 1
* @param [out] frac        basis evaluated at point \a rst_in length \a nfrac
* @param [in]  rst_in      natural coordinates input to evaluate the basis at  (-1<r,s,t<1, length 3)
* @param [in]  dimfrac     large integer (set inside Tioga to be 11^3=1331 for a p=10 hex) used for allocation frac if necessary
*/
void create_donor_frac(int *cell_index,double *xyz,int *nfrac,int *index,double *frac,double *rst_in,int *dimfrac);

void convert_to_receptor_coefficients(int *cell_index,int *npts,double *f,int *tm,int *index_out,Real *a_out);

/** Driver interface callback function:
*      callback function for Tioga that converts nodal points into solution coefficients
*      * nodal DG code this is simply a copy
*
* @param [in]  cell_index  index of cell
* @param [in]  npts        number of points to convert
* @param [in]  f           value of nodal function to be converted
* @param [out] tm          number of nodes/modes being passed out
* @param [out] index_out   index to modify solution array base 1
* @param [out] a_out       solution coefficients
void convert_to_receptor_coefficients(int *cell_index,int *npts,double *f,int *tm,int *index_out,Real *a_out);

[ ] Lastly, call the following TIOGA Interface function in order through your code:

/* preprocess and connect */
void tioga_preprocess_grids_(void);
void tioga_performconnectivity_(void);
void tioga_performconnectivity_highorder_(void);

mennodeij commented 1 year ago

@ackirby88 Wow, thank you so much for the very detailed description. I'm going to try this 👍

jsitaraman commented 1 year ago

@Andrew, thank you so much for providing the info to the user community

On Mon, Mar 6, 2023, 1:35 AM Menno Deij - van Rijswijk < @.***> wrote:

@ackirby88 https://github.com/ackirby88 Wow, thank you so much for the very detailed description. I'm going to try this 👍

— Reply to this email directly, view it on GitHub https://github.com/jsitaraman/tioga/issues/67#issuecomment-1455790320, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACH3YMGNJK4J7PZDIKPOUHDW2WVWJANCNFSM6AAAAAAUT6U6C4 . You are receiving this because you commented.Message ID: @.***>

jsitaraman / tioga

Higher-order elements: how is this supposed to work? #67