1d7b241e6Sjeremylt // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2d7b241e6Sjeremylt // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3d7b241e6Sjeremylt // reserved. See files LICENSE and NOTICE for details. 4d7b241e6Sjeremylt // 5d7b241e6Sjeremylt // This file is part of CEED, a collection of benchmarks, miniapps, software 6d7b241e6Sjeremylt // libraries and APIs for efficient high-order finite element and spectral 7d7b241e6Sjeremylt // element discretizations for exascale applications. For more information and 8d7b241e6Sjeremylt // source code availability see http://github.com/ceed. 9d7b241e6Sjeremylt // 10d7b241e6Sjeremylt // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11d7b241e6Sjeremylt // a collaborative effort of two U.S. Department of Energy organizations (Office 12d7b241e6Sjeremylt // of Science and the National Nuclear Security Administration) responsible for 13d7b241e6Sjeremylt // the planning and preparation of a capable exascale ecosystem, including 14d7b241e6Sjeremylt // software, applications, hardware, advanced system engineering and early 15d7b241e6Sjeremylt // testbed platforms, in support of the nation's exascale computing imperative. 16d7b241e6Sjeremylt 17d7b241e6Sjeremylt #include <ceed-impl.h> 18d863ab9bSjeremylt #include <ceed-backend.h> 19d7b241e6Sjeremylt 20d7b241e6Sjeremylt /// @file 21d7b241e6Sjeremylt /// Implementation of public CeedElemRestriction interfaces 22d7b241e6Sjeremylt /// 23dfdf5a53Sjeremylt /// @addtogroup CeedElemRestriction 24d7b241e6Sjeremylt /// @{ 25d7b241e6Sjeremylt 26d7b241e6Sjeremylt /** 27b11c1e72Sjeremylt @brief Create a CeedElemRestriction 28d7b241e6Sjeremylt 29b11c1e72Sjeremylt @param ceed A Ceed object where the CeedElemRestriction will be created 3061dbc9d2Sjeremylt @param imode Ordering of the ncomp components, i.e. it specifies 31a8d32208Sjeremylt the ordering of the components of the L-vector used 3261dbc9d2Sjeremylt by this CeedElemRestriction. CEED_NONINTERLACED indicates 3361dbc9d2Sjeremylt the component is the outermost index and CEED_INTERLACED 34a8d32208Sjeremylt indicates the component is the innermost index in 35a8d32208Sjeremylt ordering of the L-vector. 36b11c1e72Sjeremylt @param nelem Number of elements described in the @a indices array 37b11c1e72Sjeremylt @param elemsize Size (number of "nodes") per element 388795c945Sjeremylt @param nnodes The number of nodes in the L-vector. The input CeedVector 398795c945Sjeremylt to which the restriction will be applied is of size 408795c945Sjeremylt @a nnodes * @a ncomp. This size may include data 41d7b241e6Sjeremylt used by other CeedElemRestriction objects describing 42d7b241e6Sjeremylt different types of elements. 43b11c1e72Sjeremylt @param ncomp Number of field components per interpolation node 4495bb1877Svaleriabarra (1 for scalar fields) 45b11c1e72Sjeremylt @param mtype Memory type of the @a indices array, see CeedMemType 46b11c1e72Sjeremylt @param cmode Copy mode for the @a indices array, see CeedCopyMode 478795c945Sjeremylt @param indices Array of shape [@a nelem, @a elemsize]. Row i holds the 488795c945Sjeremylt ordered list of the indices (into the input CeedVector) 498795c945Sjeremylt for the unknowns corresponding to element i, where 5034138859Sjeremylt 0 <= i < @a nelem. All indices must be in the range 5123e8ed12Sjeremylt [0, @a nnodes - 1]. 524ce2993fSjeremylt @param[out] rstr Address of the variable where the newly created 53b11c1e72Sjeremylt CeedElemRestriction will be stored 54d7b241e6Sjeremylt 55b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 56dfdf5a53Sjeremylt 57dfdf5a53Sjeremylt @ref Basic 58b11c1e72Sjeremylt **/ 5961dbc9d2Sjeremylt int CeedElemRestrictionCreate(Ceed ceed, CeedInterlaceMode imode, 60*7509a596Sjeremylt CeedInt nelem, CeedInt elemsize, CeedInt nnodes, 61*7509a596Sjeremylt CeedInt ncomp, CeedMemType mtype, 62*7509a596Sjeremylt CeedCopyMode cmode, const CeedInt *indices, 634ce2993fSjeremylt CeedElemRestriction *rstr) { 64d7b241e6Sjeremylt int ierr; 65d7b241e6Sjeremylt 665fe0d4faSjeremylt if (!ceed->ElemRestrictionCreate) { 675fe0d4faSjeremylt Ceed delegate; 68aefd8378Sjeremylt ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 69aefd8378Sjeremylt CeedChk(ierr); 705fe0d4faSjeremylt 715fe0d4faSjeremylt if (!delegate) 72c042f62fSJeremy L Thompson // LCOV_EXCL_START 73d7b241e6Sjeremylt return CeedError(ceed, 1, "Backend does not support ElemRestrictionCreate"); 74c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 755fe0d4faSjeremylt 7661dbc9d2Sjeremylt ierr = CeedElemRestrictionCreate(delegate, imode, nelem, elemsize, 778795c945Sjeremylt nnodes, ncomp, mtype, cmode, 784ce2993fSjeremylt indices, rstr); CeedChk(ierr); 795fe0d4faSjeremylt return 0; 805fe0d4faSjeremylt } 815fe0d4faSjeremylt 824ce2993fSjeremylt ierr = CeedCalloc(1, rstr); CeedChk(ierr); 834ce2993fSjeremylt (*rstr)->ceed = ceed; 84d7b241e6Sjeremylt ceed->refcount++; 854ce2993fSjeremylt (*rstr)->refcount = 1; 8661dbc9d2Sjeremylt (*rstr)->imode = imode; 874ce2993fSjeremylt (*rstr)->nelem = nelem; 884ce2993fSjeremylt (*rstr)->elemsize = elemsize; 898795c945Sjeremylt (*rstr)->nnodes = nnodes; 904ce2993fSjeremylt (*rstr)->ncomp = ncomp; 914ce2993fSjeremylt (*rstr)->nblk = nelem; 924ce2993fSjeremylt (*rstr)->blksize = 1; 934ce2993fSjeremylt ierr = ceed->ElemRestrictionCreate(mtype, cmode, indices, *rstr); CeedChk(ierr); 94d7b241e6Sjeremylt return 0; 95d7b241e6Sjeremylt } 96d7b241e6Sjeremylt 97d7b241e6Sjeremylt /** 98*7509a596Sjeremylt @brief Create a strided CeedElemRestriction 99d7b241e6Sjeremylt 100b11c1e72Sjeremylt @param ceed A Ceed object where the CeedElemRestriction will be created 101*7509a596Sjeremylt @param nelem Number of elements described by the restriction 102b11c1e72Sjeremylt @param elemsize Size (number of "nodes") per element 1038795c945Sjeremylt @param nnodes The number of nodes in the L-vector. The input CeedVector 1048795c945Sjeremylt to which the restriction will be applied is of size 1058795c945Sjeremylt @a nnodes * @a ncomp. This size may include data 106d7b241e6Sjeremylt used by other CeedElemRestriction objects describing 1078795c945Sjeremylt different types of elements. 108b11c1e72Sjeremylt @param ncomp Number of field components per interpolation node 10995bb1877Svaleriabarra (1 for scalar fields) 110*7509a596Sjeremylt @param strides Array for strides between [nodes, components, elements]. 111*7509a596Sjeremylt The data for node i, component j, element k in the 112*7509a596Sjeremylt L-vector is given by 113*7509a596Sjeremylt i*strides[0] + j*strides[1] + k*strides[2] 1144ce2993fSjeremylt @param rstr Address of the variable where the newly created 115b11c1e72Sjeremylt CeedElemRestriction will be stored 116d7b241e6Sjeremylt 117b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 118dfdf5a53Sjeremylt 119dfdf5a53Sjeremylt @ref Basic 120b11c1e72Sjeremylt **/ 121*7509a596Sjeremylt int CeedElemRestrictionCreateStrided(Ceed ceed, CeedInt nelem, CeedInt elemsize, 122a8d32208Sjeremylt CeedInt nnodes, CeedInt ncomp, 123*7509a596Sjeremylt CeedInt strides[3], 124f90c8643Sjeremylt CeedElemRestriction *rstr) { 125d7b241e6Sjeremylt int ierr; 126d7b241e6Sjeremylt 1275fe0d4faSjeremylt if (!ceed->ElemRestrictionCreate) { 1285fe0d4faSjeremylt Ceed delegate; 129aefd8378Sjeremylt ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 130aefd8378Sjeremylt CeedChk(ierr); 1315fe0d4faSjeremylt 1325fe0d4faSjeremylt if (!delegate) 133c042f62fSJeremy L Thompson // LCOV_EXCL_START 1341d102b48SJeremy L Thompson return CeedError(ceed, 1, "Backend does not support ElemRestrictionCreate"); 135c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 1365fe0d4faSjeremylt 137*7509a596Sjeremylt ierr = CeedElemRestrictionCreateStrided(delegate, nelem, elemsize, nnodes, 138*7509a596Sjeremylt ncomp, strides, rstr); CeedChk(ierr); 1395fe0d4faSjeremylt return 0; 1405fe0d4faSjeremylt } 1415fe0d4faSjeremylt 1424ce2993fSjeremylt ierr = CeedCalloc(1, rstr); CeedChk(ierr); 1434ce2993fSjeremylt (*rstr)->ceed = ceed; 144d7b241e6Sjeremylt ceed->refcount++; 1454ce2993fSjeremylt (*rstr)->refcount = 1; 1464ce2993fSjeremylt (*rstr)->nelem = nelem; 1474ce2993fSjeremylt (*rstr)->elemsize = elemsize; 1488795c945Sjeremylt (*rstr)->nnodes = nnodes; 1494ce2993fSjeremylt (*rstr)->ncomp = ncomp; 1504ce2993fSjeremylt (*rstr)->nblk = nelem; 1514ce2993fSjeremylt (*rstr)->blksize = 1; 152*7509a596Sjeremylt ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 153*7509a596Sjeremylt for (int i = 0; i<3; i++) 154*7509a596Sjeremylt (*rstr)->strides[i] = strides[i]; 1551dfeef1dSjeremylt ierr = ceed->ElemRestrictionCreate(CEED_MEM_HOST, CEED_OWN_POINTER, NULL, 1561dfeef1dSjeremylt *rstr); 1574b8bea3bSJed Brown CeedChk(ierr); 158d7b241e6Sjeremylt return 0; 159d7b241e6Sjeremylt } 160d7b241e6Sjeremylt 161d7b241e6Sjeremylt /** 162b11c1e72Sjeremylt @brief Permute and pad indices for a blocked restriction 163d7b241e6Sjeremylt 1648795c945Sjeremylt @param indices Array of shape [@a nelem, @a elemsize]. Row i holds the 1658795c945Sjeremylt ordered list of the indices (into the input CeedVector) 1668795c945Sjeremylt for the unknowns corresponding to element i, where 16734138859Sjeremylt 0 <= i < @a nelem. All indices must be in the range 1688795c945Sjeremylt [0, @a nnodes). 169ecf6354eSJed Brown @param blkindices Array of permuted and padded indices of 170ecf6354eSJed Brown shape [@a nblk, @a elemsize, @a blksize]. 171d7b241e6Sjeremylt @param nblk Number of blocks 172d7b241e6Sjeremylt @param nelem Number of elements 173d7b241e6Sjeremylt @param blksize Number of elements in a block 174d7b241e6Sjeremylt @param elemsize Size of each element 175d7b241e6Sjeremylt 176b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 177b11c1e72Sjeremylt 178dfdf5a53Sjeremylt @ref Utility 179b11c1e72Sjeremylt **/ 180dfdf5a53Sjeremylt int CeedPermutePadIndices(const CeedInt *indices, CeedInt *blkindices, 181692c2638Sjeremylt CeedInt nblk, CeedInt nelem, CeedInt blksize, 182692c2638Sjeremylt CeedInt elemsize) { 183d7b241e6Sjeremylt for (CeedInt e = 0; e < nblk*blksize; e+=blksize) 184d7b241e6Sjeremylt for (int j = 0; j < blksize; j++) 185d7b241e6Sjeremylt for (int k = 0; k < elemsize; k++) 186d7b241e6Sjeremylt blkindices[e*elemsize + k*blksize + j] 187d7b241e6Sjeremylt = indices[CeedIntMin(e+j,nelem-1)*elemsize + k]; 188dfdf5a53Sjeremylt return 0; 189d7b241e6Sjeremylt } 190d7b241e6Sjeremylt 191d7b241e6Sjeremylt /** 192b11c1e72Sjeremylt @brief Create a blocked CeedElemRestriction, typically only called by backends 193d7b241e6Sjeremylt 194d7b241e6Sjeremylt @param ceed A Ceed object where the CeedElemRestriction will be created. 19561dbc9d2Sjeremylt @param imode Ordering of the ncomp components, i.e. it specifies 196a8d32208Sjeremylt the ordering of the components of the L-vector used 19761dbc9d2Sjeremylt by this CeedElemRestriction. CEED_NONINTERLACED indicates 19861dbc9d2Sjeremylt the component is the outermost index and CEED_INTERLACED 199a8d32208Sjeremylt indicates the component is the innermost index in 200a8d32208Sjeremylt ordering of the L-vector. 201d7b241e6Sjeremylt @param nelem Number of elements described in the @a indices array. 202b11c1e72Sjeremylt @param elemsize Size (number of unknowns) per element 203b11c1e72Sjeremylt @param blksize Number of elements in a block 2048795c945Sjeremylt @param nnodes The number of nodes in the L-vector. The input CeedVector 2058795c945Sjeremylt to which the restriction will be applied is of size 2068795c945Sjeremylt @a nnodes * @a ncomp. This size may include data 207d7b241e6Sjeremylt used by other CeedElemRestriction objects describing 208d7b241e6Sjeremylt different types of elements. 20995bb1877Svaleriabarra @param ncomp Number of field components per interpolation node 21095bb1877Svaleriabarra (1 for scalar fields) 211b11c1e72Sjeremylt @param mtype Memory type of the @a indices array, see CeedMemType 212b11c1e72Sjeremylt @param cmode Copy mode for the @a indices array, see CeedCopyMode 2138795c945Sjeremylt @param indices Array of shape [@a nelem, @a elemsize]. Row i holds the 2148795c945Sjeremylt ordered list of the indices (into the input CeedVector) 2158795c945Sjeremylt for the unknowns corresponding to element i, where 21634138859Sjeremylt 0 <= i < @a nelem. All indices must be in the range 2178795c945Sjeremylt [0, @a nnodes). The backend will permute and pad this 2188795c945Sjeremylt array to the desired ordering for the blocksize, which is 2198795c945Sjeremylt typically given by the backend. The default reordering is 2208795c945Sjeremylt to interlace elements. 2214ce2993fSjeremylt @param rstr Address of the variable where the newly created 222b11c1e72Sjeremylt CeedElemRestriction will be stored 223d7b241e6Sjeremylt 224b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 225dfdf5a53Sjeremylt 226dfdf5a53Sjeremylt @ref Advanced 227b11c1e72Sjeremylt **/ 22861dbc9d2Sjeremylt int CeedElemRestrictionCreateBlocked(Ceed ceed, CeedInterlaceMode imode, 229a8d32208Sjeremylt CeedInt nelem, CeedInt elemsize, 2308795c945Sjeremylt CeedInt blksize, CeedInt nnodes, 2318795c945Sjeremylt CeedInt ncomp, CeedMemType mtype, 2328795c945Sjeremylt CeedCopyMode cmode, const CeedInt *indices, 2334ce2993fSjeremylt CeedElemRestriction *rstr) { 234d7b241e6Sjeremylt int ierr; 235d7b241e6Sjeremylt CeedInt *blkindices; 236d7b241e6Sjeremylt CeedInt nblk = (nelem / blksize) + !!(nelem % blksize); 237d7b241e6Sjeremylt 2385fe0d4faSjeremylt if (!ceed->ElemRestrictionCreateBlocked) { 2395fe0d4faSjeremylt Ceed delegate; 240aefd8378Sjeremylt ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 241aefd8378Sjeremylt CeedChk(ierr); 2425fe0d4faSjeremylt 2435fe0d4faSjeremylt if (!delegate) 244c042f62fSJeremy L Thompson // LCOV_EXCL_START 2451d102b48SJeremy L Thompson return CeedError(ceed, 1, "Backend does not support " 2461d102b48SJeremy L Thompson "ElemRestrictionCreateBlocked"); 247c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 2485fe0d4faSjeremylt 24961dbc9d2Sjeremylt ierr = CeedElemRestrictionCreateBlocked(delegate, imode, nelem, elemsize, 2508795c945Sjeremylt blksize, nnodes, ncomp, mtype, cmode, 2514ce2993fSjeremylt indices, rstr); CeedChk(ierr); 2525fe0d4faSjeremylt return 0; 2535fe0d4faSjeremylt } 254d7b241e6Sjeremylt 2554ce2993fSjeremylt ierr = CeedCalloc(1, rstr); CeedChk(ierr); 256d7b241e6Sjeremylt 257d7b241e6Sjeremylt if (indices) { 258de686571SJeremy L Thompson ierr = CeedCalloc(nblk*blksize*elemsize, &blkindices); CeedChk(ierr); 2594b8bea3bSJed Brown ierr = CeedPermutePadIndices(indices, blkindices, nblk, nelem, blksize, 2604b8bea3bSJed Brown elemsize); 261dfdf5a53Sjeremylt CeedChk(ierr); 262d7b241e6Sjeremylt } else { 263d7b241e6Sjeremylt blkindices = NULL; 264d7b241e6Sjeremylt } 265d7b241e6Sjeremylt 2664ce2993fSjeremylt (*rstr)->ceed = ceed; 267d7b241e6Sjeremylt ceed->refcount++; 2684ce2993fSjeremylt (*rstr)->refcount = 1; 26961dbc9d2Sjeremylt (*rstr)->imode = imode; 2704ce2993fSjeremylt (*rstr)->nelem = nelem; 2714ce2993fSjeremylt (*rstr)->elemsize = elemsize; 2728795c945Sjeremylt (*rstr)->nnodes = nnodes; 2734ce2993fSjeremylt (*rstr)->ncomp = ncomp; 2744ce2993fSjeremylt (*rstr)->nblk = nblk; 2754ce2993fSjeremylt (*rstr)->blksize = blksize; 276667bc5fcSjeremylt ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 2777f823360Sjeremylt (const CeedInt *) blkindices, *rstr); CeedChk(ierr); 278d7b241e6Sjeremylt 2791d102b48SJeremy L Thompson if (cmode == CEED_OWN_POINTER) { 280d7b241e6Sjeremylt ierr = CeedFree(&indices); CeedChk(ierr); 2811d102b48SJeremy L Thompson } 282d7b241e6Sjeremylt 283d7b241e6Sjeremylt return 0; 284d7b241e6Sjeremylt } 285d7b241e6Sjeremylt 286b11c1e72Sjeremylt /** 287*7509a596Sjeremylt @brief Create a blocked strided CeedElemRestriction 288*7509a596Sjeremylt 289*7509a596Sjeremylt @param ceed A Ceed object where the CeedElemRestriction will be created 290*7509a596Sjeremylt @param nelem Number of elements described by the restriction 291*7509a596Sjeremylt @param elemsize Size (number of "nodes") per element 292*7509a596Sjeremylt @param blksize Number of elements in a block 293*7509a596Sjeremylt @param nnodes The number of nodes in the L-vector. The input CeedVector 294*7509a596Sjeremylt to which the restriction will be applied is of size 295*7509a596Sjeremylt @a nnodes * @a ncomp. This size may include data 296*7509a596Sjeremylt used by other CeedElemRestriction objects describing 297*7509a596Sjeremylt different types of elements. 298*7509a596Sjeremylt @param ncomp Number of field components per interpolation node 299*7509a596Sjeremylt (1 for scalar fields) 300*7509a596Sjeremylt @param strides Array for strides between [nodes, components, elements]. 301*7509a596Sjeremylt The data for node i, component j, element k in the 302*7509a596Sjeremylt L-vector is given by 303*7509a596Sjeremylt i*strides[0] + j*strides[1] + k*strides[2] 304*7509a596Sjeremylt @param rstr Address of the variable where the newly created 305*7509a596Sjeremylt CeedElemRestriction will be stored 306*7509a596Sjeremylt 307*7509a596Sjeremylt @return An error code: 0 - success, otherwise - failure 308*7509a596Sjeremylt 309*7509a596Sjeremylt @ref Basic 310*7509a596Sjeremylt **/ 311*7509a596Sjeremylt int CeedElemRestrictionCreateBlockedStrided(Ceed ceed, CeedInt nelem, 312*7509a596Sjeremylt CeedInt elemsize, CeedInt blksize, 313*7509a596Sjeremylt CeedInt nnodes, CeedInt ncomp, 314*7509a596Sjeremylt CeedInt strides[3], 315*7509a596Sjeremylt CeedElemRestriction *rstr) { 316*7509a596Sjeremylt int ierr; 317*7509a596Sjeremylt CeedInt nblk = (nelem / blksize) + !!(nelem % blksize); 318*7509a596Sjeremylt 319*7509a596Sjeremylt if (!ceed->ElemRestrictionCreateBlocked) { 320*7509a596Sjeremylt Ceed delegate; 321*7509a596Sjeremylt ierr = CeedGetObjectDelegate(ceed, &delegate, "ElemRestriction"); 322*7509a596Sjeremylt CeedChk(ierr); 323*7509a596Sjeremylt 324*7509a596Sjeremylt if (!delegate) 325*7509a596Sjeremylt // LCOV_EXCL_START 326*7509a596Sjeremylt return CeedError(ceed, 1, "Backend does not support " 327*7509a596Sjeremylt "ElemRestrictionCreateBlocked"); 328*7509a596Sjeremylt // LCOV_EXCL_STOP 329*7509a596Sjeremylt 330*7509a596Sjeremylt ierr = CeedElemRestrictionCreateBlockedStrided(delegate, nelem, elemsize, 331*7509a596Sjeremylt blksize, nnodes, ncomp, 332*7509a596Sjeremylt strides, rstr); 333*7509a596Sjeremylt CeedChk(ierr); 334*7509a596Sjeremylt return 0; 335*7509a596Sjeremylt } 336*7509a596Sjeremylt 337*7509a596Sjeremylt ierr = CeedCalloc(1, rstr); CeedChk(ierr); 338*7509a596Sjeremylt 339*7509a596Sjeremylt (*rstr)->ceed = ceed; 340*7509a596Sjeremylt ceed->refcount++; 341*7509a596Sjeremylt (*rstr)->refcount = 1; 342*7509a596Sjeremylt (*rstr)->nelem = nelem; 343*7509a596Sjeremylt (*rstr)->elemsize = elemsize; 344*7509a596Sjeremylt (*rstr)->nnodes = nnodes; 345*7509a596Sjeremylt (*rstr)->ncomp = ncomp; 346*7509a596Sjeremylt (*rstr)->nblk = nblk; 347*7509a596Sjeremylt (*rstr)->blksize = blksize; 348*7509a596Sjeremylt ierr = CeedMalloc(3, &(*rstr)->strides); CeedChk(ierr); 349*7509a596Sjeremylt for (int i = 0; i<3; i++) 350*7509a596Sjeremylt (*rstr)->strides[i] = strides[i]; 351*7509a596Sjeremylt ierr = ceed->ElemRestrictionCreateBlocked(CEED_MEM_HOST, CEED_OWN_POINTER, 352*7509a596Sjeremylt NULL, *rstr); CeedChk(ierr); 353*7509a596Sjeremylt 354*7509a596Sjeremylt return 0; 355*7509a596Sjeremylt } 356*7509a596Sjeremylt 357*7509a596Sjeremylt /** 358b11c1e72Sjeremylt @brief Create CeedVectors associated with a CeedElemRestriction 359b11c1e72Sjeremylt 3604ce2993fSjeremylt @param rstr CeedElemRestriction 361b11c1e72Sjeremylt @param lvec The address of the L-vector to be created, or NULL 362b11c1e72Sjeremylt @param evec The address of the E-vector to be created, or NULL 363b11c1e72Sjeremylt 364b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 365dfdf5a53Sjeremylt 366dfdf5a53Sjeremylt @ref Advanced 367b11c1e72Sjeremylt **/ 3684ce2993fSjeremylt int CeedElemRestrictionCreateVector(CeedElemRestriction rstr, CeedVector *lvec, 369d7b241e6Sjeremylt CeedVector *evec) { 370d7b241e6Sjeremylt int ierr; 371d7b241e6Sjeremylt CeedInt n, m; 3728795c945Sjeremylt m = rstr->nnodes * rstr->ncomp; 3734ce2993fSjeremylt n = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 374d7b241e6Sjeremylt if (lvec) { 3754ce2993fSjeremylt ierr = CeedVectorCreate(rstr->ceed, m, lvec); CeedChk(ierr); 376d7b241e6Sjeremylt } 377d7b241e6Sjeremylt if (evec) { 3784ce2993fSjeremylt ierr = CeedVectorCreate(rstr->ceed, n, evec); CeedChk(ierr); 379d7b241e6Sjeremylt } 380d7b241e6Sjeremylt return 0; 381d7b241e6Sjeremylt } 382d7b241e6Sjeremylt 383d7b241e6Sjeremylt /** 384d9e1f99aSValeria Barra @brief Restrict an L-vector to an E-vector or apply its transpose 385d7b241e6Sjeremylt 3864ce2993fSjeremylt @param rstr CeedElemRestriction 387d7b241e6Sjeremylt @param tmode Apply restriction or transpose 3887aaeacdcSjeremylt @param u Input vector (of shape [@a nnodes, @a ncomp] when 38961dbc9d2Sjeremylt tmode=CEED_NOTRANSPOSE, imode=CEED_INTERLACED) 390a8d32208Sjeremylt @param ru Output vector (of shape [@a nelem * @a elemsize] when 3917aaeacdcSjeremylt tmode=CEED_NOTRANSPOSE). Ordering of the e-vector is decided 3927aaeacdcSjeremylt by the backend. 393d7b241e6Sjeremylt @param request Request or CEED_REQUEST_IMMEDIATE 394b11c1e72Sjeremylt 395b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 396dfdf5a53Sjeremylt 397dfdf5a53Sjeremylt @ref Advanced 398b11c1e72Sjeremylt **/ 3994ce2993fSjeremylt int CeedElemRestrictionApply(CeedElemRestriction rstr, CeedTransposeMode tmode, 400a8d32208Sjeremylt CeedVector u, CeedVector ru, 401a8d32208Sjeremylt CeedRequest *request) { 402d7b241e6Sjeremylt CeedInt m,n; 403d7b241e6Sjeremylt int ierr; 404d7b241e6Sjeremylt 405d7b241e6Sjeremylt if (tmode == CEED_NOTRANSPOSE) { 4064ce2993fSjeremylt m = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 4078795c945Sjeremylt n = rstr->nnodes * rstr->ncomp; 408d7b241e6Sjeremylt } else { 4098795c945Sjeremylt m = rstr->nnodes * rstr->ncomp; 4104ce2993fSjeremylt n = rstr->nblk * rstr->blksize * rstr->elemsize * rstr->ncomp; 411d7b241e6Sjeremylt } 412d7b241e6Sjeremylt if (n != u->length) 413c042f62fSJeremy L Thompson // LCOV_EXCL_START 4141d102b48SJeremy L Thompson return CeedError(rstr->ceed, 2, "Input vector size %d not compatible with " 4151d102b48SJeremy L Thompson "element restriction (%d, %d)", u->length, m, n); 416c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 417a8d32208Sjeremylt if (m != ru->length) 418c042f62fSJeremy L Thompson // LCOV_EXCL_START 4191d102b48SJeremy L Thompson return CeedError(rstr->ceed, 2, "Output vector size %d not compatible with " 420a8d32208Sjeremylt "element restriction (%d, %d)", ru->length, m, n); 421c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 422074cb416Sjeremylt ierr = rstr->Apply(rstr, tmode, u, ru, request); CeedChk(ierr); 423d7b241e6Sjeremylt 424d7b241e6Sjeremylt return 0; 425d7b241e6Sjeremylt } 426d7b241e6Sjeremylt 427d7b241e6Sjeremylt /** 428d9e1f99aSValeria Barra @brief Restrict an L-vector to a block of an E-vector or apply its transpose 429be9261b7Sjeremylt 430be9261b7Sjeremylt @param rstr CeedElemRestriction 4311f37b403Sjeremylt @param block Block number to restrict to/from, i.e. block=0 will handle 4321f37b403Sjeremylt elements [0 : blksize] and block=3 will handle elements 4331f37b403Sjeremylt [3*blksize : 4*blksize] 434be9261b7Sjeremylt @param tmode Apply restriction or transpose 4357aaeacdcSjeremylt @param u Input vector (of shape [@a nnodes, @a ncomp] when 43661dbc9d2Sjeremylt tmode=CEED_NOTRANSPOSE, imode=CEED_INTERLACED) 437a8d32208Sjeremylt @param ru Output vector (of shape [@a blksize * @a elemsize] when 4387aaeacdcSjeremylt tmode=CEED_NOTRANSPOSE). Ordering of the e-vector is decided 4397aaeacdcSjeremylt by the backend. 440be9261b7Sjeremylt @param request Request or CEED_REQUEST_IMMEDIATE 441be9261b7Sjeremylt 442be9261b7Sjeremylt @return An error code: 0 - success, otherwise - failure 443be9261b7Sjeremylt 444be9261b7Sjeremylt @ref Advanced 445be9261b7Sjeremylt **/ 446be9261b7Sjeremylt int CeedElemRestrictionApplyBlock(CeedElemRestriction rstr, CeedInt block, 447a8d32208Sjeremylt CeedTransposeMode tmode, CeedVector u, 448a8d32208Sjeremylt CeedVector ru, CeedRequest *request) { 449be9261b7Sjeremylt CeedInt m,n; 450be9261b7Sjeremylt int ierr; 451be9261b7Sjeremylt 452be9261b7Sjeremylt if (tmode == CEED_NOTRANSPOSE) { 453be9261b7Sjeremylt m = rstr->blksize * rstr->elemsize * rstr->ncomp; 4548795c945Sjeremylt n = rstr->nnodes * rstr->ncomp; 455be9261b7Sjeremylt } else { 4568795c945Sjeremylt m = rstr->nnodes * rstr->ncomp; 457be9261b7Sjeremylt n = rstr->blksize * rstr->elemsize * rstr->ncomp; 458be9261b7Sjeremylt } 459be9261b7Sjeremylt if (n != u->length) 460c042f62fSJeremy L Thompson // LCOV_EXCL_START 4611d102b48SJeremy L Thompson return CeedError(rstr->ceed, 2, "Input vector size %d not compatible with " 4621d102b48SJeremy L Thompson "element restriction (%d, %d)", u->length, m, n); 463c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 464a8d32208Sjeremylt if (m != ru->length) 465c042f62fSJeremy L Thompson // LCOV_EXCL_START 4661d102b48SJeremy L Thompson return CeedError(rstr->ceed, 2, "Output vector size %d not compatible with " 467a8d32208Sjeremylt "element restriction (%d, %d)", ru->length, m, n); 468c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 469be9261b7Sjeremylt if (rstr->blksize*block > rstr->nelem) 470c042f62fSJeremy L Thompson // LCOV_EXCL_START 4711d102b48SJeremy L Thompson return CeedError(rstr->ceed, 2, "Cannot retrieve block %d, element %d > " 4721d102b48SJeremy L Thompson "total elements %d", block, rstr->blksize*block, 4731d102b48SJeremy L Thompson rstr->nelem); 474c042f62fSJeremy L Thompson // LCOV_EXCL_STOP 475074cb416Sjeremylt ierr = rstr->ApplyBlock(rstr, block, tmode, u, ru, request); 476be9261b7Sjeremylt CeedChk(ierr); 477be9261b7Sjeremylt 478be9261b7Sjeremylt return 0; 479be9261b7Sjeremylt } 480be9261b7Sjeremylt 481be9261b7Sjeremylt /** 482d9e1f99aSValeria Barra @brief Get the multiplicity of nodes in a CeedElemRestriction 4831469ee4dSjeremylt 4841469ee4dSjeremylt @param rstr CeedElemRestriction 485b9c05c73SJed Brown @param[out] mult Vector to store multiplicity (of size nnodes*ncomp) 4861469ee4dSjeremylt 4871469ee4dSjeremylt @return An error code: 0 - success, otherwise - failure 4881469ee4dSjeremylt 4891469ee4dSjeremylt @ref Advanced 4901469ee4dSjeremylt **/ 4911469ee4dSjeremylt int CeedElemRestrictionGetMultiplicity(CeedElemRestriction rstr, 4921469ee4dSjeremylt CeedVector mult) { 4931469ee4dSjeremylt int ierr; 4941469ee4dSjeremylt CeedVector evec; 4951469ee4dSjeremylt 4961469ee4dSjeremylt // Create and set evec 4971469ee4dSjeremylt ierr = CeedElemRestrictionCreateVector(rstr, NULL, &evec); CeedChk(ierr); 4981469ee4dSjeremylt ierr = CeedVectorSetValue(evec, 1.0); CeedChk(ierr); 499fa9eac48SJed Brown ierr = CeedVectorSetValue(mult, 0.0); CeedChk(ierr); 5001469ee4dSjeremylt 5011469ee4dSjeremylt // Apply to get multiplicity 502a8d32208Sjeremylt ierr = CeedElemRestrictionApply(rstr, CEED_TRANSPOSE, evec, mult, 503efc78312Sjeremylt CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 5041469ee4dSjeremylt 5051469ee4dSjeremylt // Cleanup 5061469ee4dSjeremylt ierr = CeedVectorDestroy(&evec); CeedChk(ierr); 5071469ee4dSjeremylt 5081469ee4dSjeremylt return 0; 5091469ee4dSjeremylt } 5101469ee4dSjeremylt 5111469ee4dSjeremylt /** 5124ce2993fSjeremylt @brief Get the Ceed associated with a CeedElemRestriction 5134ce2993fSjeremylt 5144ce2993fSjeremylt @param rstr CeedElemRestriction 5154ce2993fSjeremylt @param[out] ceed Variable to store Ceed 5164ce2993fSjeremylt 5174ce2993fSjeremylt @return An error code: 0 - success, otherwise - failure 5184ce2993fSjeremylt 51923617272Sjeremylt @ref Advanced 5204ce2993fSjeremylt **/ 5214ce2993fSjeremylt int CeedElemRestrictionGetCeed(CeedElemRestriction rstr, Ceed *ceed) { 5224ce2993fSjeremylt *ceed = rstr->ceed; 5234ce2993fSjeremylt return 0; 5244ce2993fSjeremylt } 5254ce2993fSjeremylt 5264ce2993fSjeremylt /** 52761dbc9d2Sjeremylt @brief Get the L-vector interlaced mode of a CeedElemRestriction 528a8d32208Sjeremylt 529a8d32208Sjeremylt @param rstr CeedElemRestriction 53061dbc9d2Sjeremylt @param[out] imode Variable to store imode 531a8d32208Sjeremylt 532a8d32208Sjeremylt @return An error code: 0 - success, otherwise - failure 533a8d32208Sjeremylt 534a8d32208Sjeremylt @ref Advanced 535a8d32208Sjeremylt **/ 53661dbc9d2Sjeremylt int CeedElemRestrictionGetIMode(CeedElemRestriction rstr, 53761dbc9d2Sjeremylt CeedInterlaceMode *imode) { 538*7509a596Sjeremylt if (rstr->strides) 539*7509a596Sjeremylt // LCOV_EXCL_START 540*7509a596Sjeremylt return CeedError(rstr->ceed, 1, "Strided ElemRestriction has no interlace " 541*7509a596Sjeremylt "mode"); 542*7509a596Sjeremylt // LCOV_EXCL_STOP 543*7509a596Sjeremylt 54461dbc9d2Sjeremylt *imode = rstr->imode; 545a8d32208Sjeremylt return 0; 546a8d32208Sjeremylt } 547a8d32208Sjeremylt 548a8d32208Sjeremylt /** 549b11c1e72Sjeremylt @brief Get the total number of elements in the range of a CeedElemRestriction 550d7b241e6Sjeremylt 5514ce2993fSjeremylt @param rstr CeedElemRestriction 552288c0443SJeremy L Thompson @param[out] numelem Variable to store number of elements 553b11c1e72Sjeremylt 554b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 555dfdf5a53Sjeremylt 55623617272Sjeremylt @ref Advanced 557b11c1e72Sjeremylt **/ 5584ce2993fSjeremylt int CeedElemRestrictionGetNumElements(CeedElemRestriction rstr, 5594ce2993fSjeremylt CeedInt *numelem) { 5604ce2993fSjeremylt *numelem = rstr->nelem; 5614ce2993fSjeremylt return 0; 5624ce2993fSjeremylt } 5634ce2993fSjeremylt 5644ce2993fSjeremylt /** 5654ce2993fSjeremylt @brief Get the size of elements in the CeedElemRestriction 5664ce2993fSjeremylt 5674ce2993fSjeremylt @param rstr CeedElemRestriction 5684ce2993fSjeremylt @param[out] elemsize Variable to store size of elements 5694ce2993fSjeremylt 5704ce2993fSjeremylt @return An error code: 0 - success, otherwise - failure 5714ce2993fSjeremylt 57223617272Sjeremylt @ref Advanced 5734ce2993fSjeremylt **/ 5744ce2993fSjeremylt int CeedElemRestrictionGetElementSize(CeedElemRestriction rstr, 5754ce2993fSjeremylt CeedInt *elemsize) { 5764ce2993fSjeremylt *elemsize = rstr->elemsize; 5774ce2993fSjeremylt return 0; 5784ce2993fSjeremylt } 5794ce2993fSjeremylt 5804ce2993fSjeremylt /** 5814ce2993fSjeremylt @brief Get the number of degrees of freedom in the range of a 5824ce2993fSjeremylt CeedElemRestriction 5834ce2993fSjeremylt 5844ce2993fSjeremylt @param rstr CeedElemRestriction 5858795c945Sjeremylt @param[out] numnodes Variable to store number of nodes 5864ce2993fSjeremylt 5874ce2993fSjeremylt @return An error code: 0 - success, otherwise - failure 5884ce2993fSjeremylt 58923617272Sjeremylt @ref Advanced 5904ce2993fSjeremylt **/ 5918795c945Sjeremylt int CeedElemRestrictionGetNumNodes(CeedElemRestriction rstr, 5928795c945Sjeremylt CeedInt *numnodes) { 5938795c945Sjeremylt *numnodes = rstr->nnodes; 5944ce2993fSjeremylt return 0; 5954ce2993fSjeremylt } 5964ce2993fSjeremylt 5974ce2993fSjeremylt /** 5984ce2993fSjeremylt @brief Get the number of components in the elements of a 5994ce2993fSjeremylt CeedElemRestriction 6004ce2993fSjeremylt 6014ce2993fSjeremylt @param rstr CeedElemRestriction 6024ce2993fSjeremylt @param[out] numcomp Variable to store number of components 6034ce2993fSjeremylt 6044ce2993fSjeremylt @return An error code: 0 - success, otherwise - failure 6054ce2993fSjeremylt 60623617272Sjeremylt @ref Advanced 6074ce2993fSjeremylt **/ 6084ce2993fSjeremylt int CeedElemRestrictionGetNumComponents(CeedElemRestriction rstr, 6094ce2993fSjeremylt CeedInt *numcomp) { 6104ce2993fSjeremylt *numcomp = rstr->ncomp; 6114ce2993fSjeremylt return 0; 6124ce2993fSjeremylt } 6134ce2993fSjeremylt 6144ce2993fSjeremylt /** 6154ce2993fSjeremylt @brief Get the number of blocks in a CeedElemRestriction 6164ce2993fSjeremylt 6174ce2993fSjeremylt @param rstr CeedElemRestriction 6184ce2993fSjeremylt @param[out] numblock Variable to store number of blocks 6194ce2993fSjeremylt 6204ce2993fSjeremylt @return An error code: 0 - success, otherwise - failure 6214ce2993fSjeremylt 62223617272Sjeremylt @ref Advanced 6234ce2993fSjeremylt **/ 6244ce2993fSjeremylt int CeedElemRestrictionGetNumBlocks(CeedElemRestriction rstr, 6254ce2993fSjeremylt CeedInt *numblock) { 6264ce2993fSjeremylt *numblock = rstr->nblk; 6274ce2993fSjeremylt return 0; 6284ce2993fSjeremylt } 6294ce2993fSjeremylt 6304ce2993fSjeremylt /** 6314ce2993fSjeremylt @brief Get the size of blocks in the CeedElemRestriction 6324ce2993fSjeremylt 633288c0443SJeremy L Thompson @param rstr CeedElemRestriction 6344ce2993fSjeremylt @param[out] blksize Variable to store size of blocks 6354ce2993fSjeremylt 6364ce2993fSjeremylt @return An error code: 0 - success, otherwise - failure 6374ce2993fSjeremylt 63823617272Sjeremylt @ref Advanced 6394ce2993fSjeremylt **/ 6404ce2993fSjeremylt int CeedElemRestrictionGetBlockSize(CeedElemRestriction rstr, 6414ce2993fSjeremylt CeedInt *blksize) { 6424ce2993fSjeremylt *blksize = rstr->blksize; 6434ce2993fSjeremylt return 0; 6444ce2993fSjeremylt } 6454ce2993fSjeremylt 6464ce2993fSjeremylt /** 647*7509a596Sjeremylt @brief Get the strides of a strided CeedElemRestriction 648*7509a596Sjeremylt 649*7509a596Sjeremylt @param rstr CeedElemRestriction 650*7509a596Sjeremylt @param[out] strides Variable to store strides array 651*7509a596Sjeremylt 652*7509a596Sjeremylt @return An error code: 0 - success, otherwise - failure 653*7509a596Sjeremylt 654*7509a596Sjeremylt @ref Advanced 655*7509a596Sjeremylt **/ 656*7509a596Sjeremylt int CeedElemRestrictionGetStrides(CeedElemRestriction rstr, 657*7509a596Sjeremylt CeedInt (*strides)[3]) { 658*7509a596Sjeremylt if (!rstr->strides) 659*7509a596Sjeremylt // LCOV_EXCL_START 660*7509a596Sjeremylt return CeedError(rstr->ceed, 1, "ElemRestriction has no stride data"); 661*7509a596Sjeremylt // LCOV_EXCL_STOP 662*7509a596Sjeremylt 663*7509a596Sjeremylt for (int i = 0; i<3; i++) 664*7509a596Sjeremylt (*strides)[i] = rstr->strides[i]; 665*7509a596Sjeremylt return 0; 666*7509a596Sjeremylt } 667*7509a596Sjeremylt 668*7509a596Sjeremylt /** 6694ce2993fSjeremylt @brief Get the backend data of a CeedElemRestriction 6704ce2993fSjeremylt 671288c0443SJeremy L Thompson @param rstr CeedElemRestriction 6724ce2993fSjeremylt @param[out] data Variable to store data 6734ce2993fSjeremylt 6744ce2993fSjeremylt @return An error code: 0 - success, otherwise - failure 6754ce2993fSjeremylt 67623617272Sjeremylt @ref Advanced 6774ce2993fSjeremylt **/ 6781d102b48SJeremy L Thompson int CeedElemRestrictionGetData(CeedElemRestriction rstr, void **data) { 6794ce2993fSjeremylt *data = rstr->data; 680d7b241e6Sjeremylt return 0; 681d7b241e6Sjeremylt } 682d7b241e6Sjeremylt 683d7b241e6Sjeremylt /** 684fe2413ffSjeremylt @brief Set the backend data of a CeedElemRestriction 685fe2413ffSjeremylt 686288c0443SJeremy L Thompson @param[out] rstr CeedElemRestriction 687fe2413ffSjeremylt @param data Data to set 688fe2413ffSjeremylt 689fe2413ffSjeremylt @return An error code: 0 - success, otherwise - failure 690fe2413ffSjeremylt 691fe2413ffSjeremylt @ref Advanced 692fe2413ffSjeremylt **/ 6931d102b48SJeremy L Thompson int CeedElemRestrictionSetData(CeedElemRestriction rstr, void **data) { 694fe2413ffSjeremylt rstr->data = *data; 695fe2413ffSjeremylt return 0; 696fe2413ffSjeremylt } 697fe2413ffSjeremylt 698fe2413ffSjeremylt /** 699f02ca4a2SJed Brown @brief View a CeedElemRestriction 700f02ca4a2SJed Brown 701f02ca4a2SJed Brown @param[in] rstr CeedElemRestriction to view 702f02ca4a2SJed Brown @param[in] stream Stream to write; typically stdout/stderr or a file 703f02ca4a2SJed Brown 704f02ca4a2SJed Brown @return Error code: 0 - success, otherwise - failure 705f02ca4a2SJed Brown 706f02ca4a2SJed Brown @ref Utility 707f02ca4a2SJed Brown **/ 708f02ca4a2SJed Brown int CeedElemRestrictionView(CeedElemRestriction rstr, FILE *stream) { 709*7509a596Sjeremylt char stridesstr[500]; 710*7509a596Sjeremylt if (rstr->strides) 711*7509a596Sjeremylt sprintf(stridesstr, "[%d, %d, %d]", rstr->strides[0], rstr->strides[1], 712*7509a596Sjeremylt rstr->strides[2]); 713*7509a596Sjeremylt 7141d102b48SJeremy L Thompson fprintf(stream, "CeedElemRestriction from (%d, %d) to %d elements with %d " 715*7509a596Sjeremylt "nodes each and %s %s\n", rstr->nnodes, rstr->ncomp, rstr->nelem, 716*7509a596Sjeremylt rstr->elemsize, rstr->strides ? "strides" : "L-vector components", 717*7509a596Sjeremylt rstr->strides ? stridesstr : CeedInterlaceModes[rstr->imode]); 718f02ca4a2SJed Brown return 0; 719f02ca4a2SJed Brown } 720f02ca4a2SJed Brown 721f02ca4a2SJed Brown /** 722b11c1e72Sjeremylt @brief Destroy a CeedElemRestriction 723b11c1e72Sjeremylt 7244ce2993fSjeremylt @param rstr CeedElemRestriction to destroy 725b11c1e72Sjeremylt 726b11c1e72Sjeremylt @return An error code: 0 - success, otherwise - failure 727dfdf5a53Sjeremylt 728dfdf5a53Sjeremylt @ref Basic 729b11c1e72Sjeremylt **/ 7304ce2993fSjeremylt int CeedElemRestrictionDestroy(CeedElemRestriction *rstr) { 731d7b241e6Sjeremylt int ierr; 732d7b241e6Sjeremylt 7331d102b48SJeremy L Thompson if (!*rstr || --(*rstr)->refcount > 0) 7341d102b48SJeremy L Thompson return 0; 7354ce2993fSjeremylt if ((*rstr)->Destroy) { 7364ce2993fSjeremylt ierr = (*rstr)->Destroy(*rstr); CeedChk(ierr); 737d7b241e6Sjeremylt } 738*7509a596Sjeremylt ierr = CeedFree(&(*rstr)->strides); CeedChk(ierr); 7394ce2993fSjeremylt ierr = CeedDestroy(&(*rstr)->ceed); CeedChk(ierr); 7404ce2993fSjeremylt ierr = CeedFree(rstr); CeedChk(ierr); 741d7b241e6Sjeremylt return 0; 742d7b241e6Sjeremylt } 743d7b241e6Sjeremylt 744d7b241e6Sjeremylt /// @} 745