1 // Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include "../qfunctions/inverse_multiplicity.h" 9 #include "../navierstokes.h" 10 11 /** 12 * @brief Get the inverse of the node multiplicity 13 * 14 * Local multiplicity concerns only the number of elements associated with a node in the current rank's partition. 15 * Global multiplicity is the multiplicity for the global mesh, including periodicity. 16 * 17 * @param[in] ceed `Ceed` object for the output `CeedVector` and `CeedElemRestriction` 18 * @param[in] dm `DM` for the grid 19 * @param[in] domain_label `DMLabel` for `DMPlex` domain 20 * @param[in] label_value Stratum value 21 * @param[in] height Height of `DMPlex` topology 22 * @param[in] dm_field Index of `DMPlex` field 23 * @param[in] get_global_multiplicity Whether the multiplicity should be global or local 24 * @param[out] elem_restr_inv_multiplicity `CeedElemRestriction` needed to use the multiplicity vector 25 * @param[out] inv_multiplicity `CeedVector` containing the inverse of the multiplicity 26 */ 27 PetscErrorCode GetInverseMultiplicity(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt dm_field, 28 PetscBool get_global_multiplicity, CeedElemRestriction *elem_restr_inv_multiplicity, 29 CeedVector *inv_multiplicity) { 30 Vec Multiplicity, Multiplicity_loc; 31 PetscMemType m_mem_type; 32 CeedVector multiplicity; 33 CeedQFunction qf_multiplicity; 34 CeedOperator op_multiplicity; 35 CeedInt num_comp; 36 CeedElemRestriction elem_restr; 37 38 PetscFunctionBeginUser; 39 40 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, label_value, height, 1, elem_restr_inv_multiplicity)); 41 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_inv_multiplicity, inv_multiplicity, NULL)); 42 43 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, label_value, height, dm_field, &elem_restr)); 44 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr, &num_comp)); 45 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr, &multiplicity, NULL)); 46 47 if (get_global_multiplicity) { 48 // In order to get global multiplicity, we need to run through DMLocalToGlobal -> DMGlobalToLocal. 49 PetscCall(DMGetLocalVector(dm, &Multiplicity_loc)); 50 PetscCall(DMGetGlobalVector(dm, &Multiplicity)); 51 52 PetscCall(VecPetscToCeed(Multiplicity_loc, &m_mem_type, multiplicity)); 53 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr, multiplicity)); 54 PetscCall(VecCeedToPetsc(multiplicity, m_mem_type, Multiplicity_loc)); 55 PetscCall(VecZeroEntries(Multiplicity)); 56 PetscCall(DMLocalToGlobal(dm, Multiplicity_loc, ADD_VALUES, Multiplicity)); 57 PetscCall(DMGlobalToLocal(dm, Multiplicity, INSERT_VALUES, Multiplicity_loc)); 58 PetscCall(VecPetscToCeed(Multiplicity_loc, &m_mem_type, multiplicity)); 59 } else { 60 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr, multiplicity)); 61 } 62 63 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, InverseMultiplicity, InverseMultiplicity_loc, &qf_multiplicity)); 64 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_multiplicity, "multiplicity", num_comp, CEED_EVAL_NONE)); 65 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_multiplicity, "inverse multiplicity", 1, CEED_EVAL_NONE)); 66 67 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_multiplicity, NULL, NULL, &op_multiplicity)); 68 PetscCallCeed(ceed, CeedOperatorSetName(op_multiplicity, "InverseMultiplicity")); 69 PetscCallCeed(ceed, CeedOperatorSetField(op_multiplicity, "multiplicity", elem_restr, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 70 PetscCallCeed(ceed, 71 CeedOperatorSetField(op_multiplicity, "inverse multiplicity", *elem_restr_inv_multiplicity, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 72 73 PetscCallCeed(ceed, CeedOperatorApply(op_multiplicity, multiplicity, *inv_multiplicity, CEED_REQUEST_IMMEDIATE)); 74 75 if (get_global_multiplicity) { 76 PetscCall(VecCeedToPetsc(multiplicity, m_mem_type, Multiplicity_loc)); 77 PetscCall(DMRestoreLocalVector(dm, &Multiplicity_loc)); 78 PetscCall(DMRestoreGlobalVector(dm, &Multiplicity)); 79 } 80 PetscCallCeed(ceed, CeedVectorDestroy(&multiplicity)); 81 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr)); 82 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_multiplicity)); 83 PetscCallCeed(ceed, CeedOperatorDestroy(&op_multiplicity)); 84 PetscFunctionReturn(PETSC_SUCCESS); 85 } 86