1 // Copyright (c) 2017-2026, 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 */
GetInverseMultiplicity(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,PetscInt dm_field,PetscBool get_global_multiplicity,CeedElemRestriction * elem_restr_inv_multiplicity,CeedVector * inv_multiplicity)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