From f037b23fb180c8c5ce0ae4c93fe234dc6debb6e2 Mon Sep 17 00:00:00 2001 From: "Mittal, Ketan" Date: Wed, 24 Apr 2024 12:10:08 -0700 Subject: [PATCH 1/9] initial commit --- fem/gridfunc.cpp | 9 +- fem/gridfunc.hpp | 4 +- fem/pgridfunc.cpp | 2 +- fem/pgridfunc.hpp | 2 +- fem/tmop.cpp | 322 +++++++++++++++++++---------- fem/tmop.hpp | 11 +- fem/tmop_tools.cpp | 2 + miniapps/meshing/pmesh-fitting.cpp | 33 ++- 8 files changed, 262 insertions(+), 123 deletions(-) diff --git a/fem/gridfunc.cpp b/fem/gridfunc.cpp index 6b2281c88f..e9393da1a3 100644 --- a/fem/gridfunc.cpp +++ b/fem/gridfunc.cpp @@ -1321,9 +1321,9 @@ void GridFunction::ProjectVectorFieldOn(GridFunction &vec_field, int comp) } } -void GridFunction::AccumulateAndCountDerivativeValues(int comp, int der_comp, - GridFunction &der, - Array &zones_per_dof) +void GridFunction::AccumulateAndCountDerivativeValues( + int comp, int der_comp, GridFunction &der, + Array &zones_per_dof) const { FiniteElementSpace * der_fes = der.FESpace(); ElementTransformation * transf; @@ -1374,7 +1374,8 @@ void GridFunction::AccumulateAndCountDerivativeValues(int comp, int der_comp, } } -void GridFunction::GetDerivative(int comp, int der_comp, GridFunction &der) +void GridFunction::GetDerivative(int comp, int der_comp, + GridFunction &der) const { Array overlap; AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap); diff --git a/fem/gridfunc.hpp b/fem/gridfunc.hpp index 4a9c40fac9..df1ef3b687 100644 --- a/fem/gridfunc.hpp +++ b/fem/gridfunc.hpp @@ -321,7 +321,7 @@ public: @param[out] der The resulting derivative (scalar function). The FiniteElementSpace of this function must be set before the call. */ - void GetDerivative(int comp, int der_comp, GridFunction &der); + void GetDerivative(int comp, int der_comp, GridFunction &der) const; real_t GetDivergence(ElementTransformation &tr) const; @@ -443,7 +443,7 @@ protected: GetDerivative() method; see its documentation. */ void AccumulateAndCountDerivativeValues(int comp, int der_comp, GridFunction &der, - Array &zones_per_dof); + Array &zones_per_dof) const; void AccumulateAndCountBdrValues(Coefficient *coeff[], VectorCoefficient *vcoeff, diff --git a/fem/pgridfunc.cpp b/fem/pgridfunc.cpp index 0f042864b4..940535ec76 100644 --- a/fem/pgridfunc.cpp +++ b/fem/pgridfunc.cpp @@ -518,7 +518,7 @@ void ParGridFunction::CountElementsPerVDof(Array &elem_per_vdof) const } void ParGridFunction::GetDerivative(int comp, int der_comp, - ParGridFunction &der) + ParGridFunction &der) const { Array overlap; AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap); diff --git a/fem/pgridfunc.hpp b/fem/pgridfunc.hpp index 833615ebe9..63d40fef3a 100644 --- a/fem/pgridfunc.hpp +++ b/fem/pgridfunc.hpp @@ -231,7 +231,7 @@ public: void CountElementsPerVDof(Array &elem_per_vdof) const override; /// Parallel version of GridFunction::GetDerivative(); see its documentation. - void GetDerivative(int comp, int der_comp, ParGridFunction &der); + void GetDerivative(int comp, int der_comp, ParGridFunction &der) const; /** Sets the output vector @a dof_vals to the values of the degrees of freedom of element @a el. If @a el is greater than or equal to the number diff --git a/fem/tmop.cpp b/fem/tmop.cpp index a3400cd187..03df678f49 100644 --- a/fem/tmop.cpp +++ b/fem/tmop.cpp @@ -2987,12 +2987,19 @@ void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &pos, void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0, const Array &smarker, Coefficient &coeff, - AdaptivityEvaluator &ae) + AdaptivityEvaluator &ae, + AdaptivityEvaluator *aegrad, + AdaptivityEvaluator *aehess) { // To have both we must duplicate the markers. MFEM_VERIFY(surf_fit_pos == NULL, "Using both fitting approaches is not supported."); + const int dim = s0.FESpace()->GetMesh()->Dimension(); + ParMesh *pmesh = s0.ParFESpace()->GetParMesh(); + MFEM_VERIFY(pmesh->GetNodes()->Size() == dim*s0.Size(), + "Mesh and level-set polynomial order must be the same."); + delete surf_fit_gf; surf_fit_gf = new GridFunction(s0); s0.CountElementsPerVDof(surf_fit_dof_count); @@ -3000,11 +3007,80 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0, surf_fit_coeff = &coeff; surf_fit_eval = &ae; - surf_fit_eval->SetParMetaInfo(*s0.ParFESpace()->GetParMesh(), - *s0.ParFESpace()); + surf_fit_eval->SetParMetaInfo(*pmesh, *s0.ParFESpace()); surf_fit_eval->SetInitialField (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf); surf_fit_gf_bg = false; + + if (!aegrad) { return; } + + // Define spaces for grad and hessian + delete surf_fit_grad; + + ParFiniteElementSpace *fes = s0.ParFESpace(); + + const H1_FECollection *fec = dynamic_cast + (surf_fit_gf->FESpace()->FEColl()); + if (!fec) { return; } + H1_FECollection *fec_grad = new H1_FECollection(fec->GetOrder(), dim, + fec->GetBasisType()); + ParFiniteElementSpace *fes_grad = new ParFiniteElementSpace(pmesh, fec_grad, + dim); + surf_fit_grad = new GridFunction(fes_grad); + surf_fit_grad->MakeOwner(fec_grad); + + for (int d = 0; d < dim; d++) + { + ParGridFunction surf_fit_grad_comp(fes, surf_fit_grad->GetData()+d*s0.Size()); + s0.GetDerivative(1, d, surf_fit_grad_comp); + } + + surf_fit_eval_bg_grad = aegrad; + surf_fit_eval_bg_grad->SetParMetaInfo(*pmesh, *fes_grad); + surf_fit_eval_bg_grad->SetInitialField(*pmesh->GetNodes(), *surf_fit_grad); + + MFEM_VERIFY(aehess,"Specify an adaptivity evaluator for the Hessian terms" + "as well."); + delete surf_fit_hess; + H1_FECollection *fec_hess = new H1_FECollection(fec->GetOrder(), dim, + fec->GetBasisType()); + ParFiniteElementSpace *fes_hess = new ParFiniteElementSpace(pmesh, fec_hess, + dim*dim); + surf_fit_hess = new GridFunction(fes_hess); + surf_fit_hess->MakeOwner(fec_hess); + + int id = 0; + for (int d = 0; d < dim; d++) + { + for (int idir = 0; idir < dim; idir++) + { + ParGridFunction surf_fit_grad_comp(fes, + surf_fit_grad->GetData()+d*s0.Size()); + ParGridFunction surf_fit_hess_comp(fes, + surf_fit_hess->GetData()+id*s0.Size()); + surf_fit_grad_comp.GetDerivative(1, idir, surf_fit_hess_comp); + id++; + } + } + + surf_fit_eval_bg_hess = aehess; + surf_fit_eval_bg_hess->SetParMetaInfo(*pmesh, *fes_hess); + surf_fit_eval_bg_hess->SetInitialField(*pmesh->GetNodes(), *surf_fit_hess); + + surf_fit_gf_bg = true; + // Store DOF indices that are marked for fitting. Used to reduce work for + // transferring information between source/background and current mesh. + surf_fit_marker_dof_index.SetSize(0); + for (int i = 0; i < surf_fit_marker->Size(); i++) + { + if ((*surf_fit_marker)[i] == true) + { + surf_fit_marker_dof_index.Append(i); + } + } + + *surf_fit_grad = 0.0; + *surf_fit_hess = 0.0; } void TMOP_Integrator::EnableSurfaceFittingFromSource( @@ -3022,7 +3098,6 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource( // Setup for level set function delete surf_fit_gf; surf_fit_gf = new GridFunction(s0); - *surf_fit_gf = 0.0; surf_fit_marker = &smarker; surf_fit_coeff = &coeff; surf_fit_eval = &ae; @@ -4245,6 +4320,15 @@ real_t TMOP_Integrator::GetSurfaceFittingWeight() return 0.0; } +void TMOP_Integrator::GetSurfaceFittingLevelSet(GridFunction &s0) +{ + + MFEM_VERIFY(s0.Size() == surf_fit_gf->Size(), + "Provided function and internal level-set function should" + "be of same size."); + s0 = *surf_fit_gf; +} + void TMOP_Integrator::EnableNormalization(const GridFunction &x) { ComputeNormalizationEnergies(x, metric_normal, lim_normal, surf_fit_normal); @@ -4376,6 +4460,132 @@ void TMOP_Integrator::ComputeMinJac(const Vector &x, dx = detv_avg_min / dxscale; } +void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, + int new_x_ordering) +{ + if (!surf_fit_gf) { return; } + if (surf_fit_gf_bg && + dynamic_cast(surf_fit_eval) && + dynamic_cast(surf_fit_eval_bg_grad) && + dynamic_cast(surf_fit_eval_bg_hess)) + { + // Interpolate information for only DOFs marked for fitting. + const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension(); + const int cnt = surf_fit_marker_dof_index.Size(); + const int total_cnt = new_x.Size()/dim; + Vector new_x_sorted(cnt*dim); + if (new_x_ordering == 0) + { + for (int d = 0; d < dim; d++) + { + for (int i = 0; i < cnt; i++) + { + int dof_index = surf_fit_marker_dof_index[i]; + new_x_sorted(i + d*cnt) = new_x(dof_index + d*total_cnt); + } + } + } + else + { + for (int i = 0; i < cnt; i++) + { + int dof_index = surf_fit_marker_dof_index[i]; + for (int d = 0; d < dim; d++) + { + new_x_sorted(d + i*dim) = new_x(d + dof_index*dim); + } + } + } + + + Vector surf_fit_gf_int, surf_fit_grad_int, surf_fit_hess_int; + surf_fit_eval->ComputeAtNewPosition(new_x_sorted, surf_fit_gf_int, + new_x_ordering); + for (int i = 0; i < cnt; i++) + { + int dof_index = surf_fit_marker_dof_index[i]; + (*surf_fit_gf)[dof_index] = surf_fit_gf_int(i); + } + + + surf_fit_eval_bg_grad->ComputeAtNewPosition(new_x_sorted, surf_fit_grad_int, + new_x_ordering); + + // Assumes surf_fit_grad and surf_fit_gf share the same space + const int grad_dim = surf_fit_grad->VectorDim(); + const int grad_cnt = surf_fit_grad->Size()/grad_dim; + if (surf_fit_grad->FESpace()->GetOrdering() == Ordering::byNODES) + { + for (int d = 0; d < grad_dim; d++) + { + for (int i = 0; i < cnt; i++) + { + int dof_index = surf_fit_marker_dof_index[i]; + (*surf_fit_grad)[dof_index + d*grad_cnt] = + surf_fit_grad_int(i + d*cnt); + } + } + } + else + { + for (int i = 0; i < cnt; i++) + { + int dof_index = surf_fit_marker_dof_index[i]; + for (int d = 0; d < grad_dim; d++) + { + (*surf_fit_grad)[dof_index*grad_dim + d] = + surf_fit_grad_int(i*grad_dim + d); + } + } + } + + surf_fit_eval_bg_hess->ComputeAtNewPosition(new_x_sorted, surf_fit_hess_int, + new_x_ordering); + // Assumes surf_fit_hess and surf_fit_gf share the same space + const int hess_dim = surf_fit_hess->VectorDim(); + const int hess_cnt = surf_fit_hess->Size()/hess_dim; + if (surf_fit_hess->FESpace()->GetOrdering() == Ordering::byNODES) + { + for (int d = 0; d < hess_dim; d++) + { + for (int i = 0; i < cnt; i++) + { + int dof_index = surf_fit_marker_dof_index[i]; + (*surf_fit_hess)[dof_index + d*hess_cnt] = + surf_fit_hess_int(i + d*cnt); + } + } + } + else + { + for (int i = 0; i < cnt; i++) + { + int dof_index = surf_fit_marker_dof_index[i]; + for (int d = 0; d < hess_dim; d++) + { + (*surf_fit_hess)[dof_index*hess_dim + d] = + surf_fit_hess_int(i*hess_dim + d); + } + } + } + + } + else + { + surf_fit_eval->ComputeAtNewPosition(new_x, *surf_fit_gf, new_x_ordering); + if (surf_fit_eval_bg_grad) + { + surf_fit_eval_bg_grad->ComputeAtNewPosition(new_x, *surf_fit_grad, + new_x_ordering); + } + if (surf_fit_eval_bg_hess) + { + surf_fit_eval_bg_hess->ComputeAtNewPosition(new_x, *surf_fit_hess, + new_x_ordering); + } + } +} + void TMOP_Integrator:: UpdateAfterMeshPositionChange(const Vector &x_new, const FiniteElementSpace &x_fes) @@ -4409,109 +4619,7 @@ UpdateAfterMeshPositionChange(const Vector &x_new, // Update surf_fit_gf if surface fitting is enabled. if (surf_fit_gf) { - if (surf_fit_gf_bg) - { - // Interpolate information for only DOFs marked for fitting. - const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension(); - const int cnt = surf_fit_marker_dof_index.Size(); - const int total_cnt = x_new.Size()/dim; - Vector new_x_sorted(cnt*dim); - if (ordering == 0) - { - for (int d = 0; d < dim; d++) - { - for (int i = 0; i < cnt; i++) - { - int dof_index = surf_fit_marker_dof_index[i]; - new_x_sorted(i + d*cnt) = x_new(dof_index + d*total_cnt); - } - } - } - else - { - for (int i = 0; i < cnt; i++) - { - int dof_index = surf_fit_marker_dof_index[i]; - for (int d = 0; d < dim; d++) - { - new_x_sorted(d + i*dim) = x_new(d + dof_index*dim); - } - } - } - - Vector surf_fit_gf_int, surf_fit_grad_int, surf_fit_hess_int; - surf_fit_eval->ComputeAtNewPosition( - new_x_sorted, surf_fit_gf_int, ordering); - for (int i = 0; i < cnt; i++) - { - int dof_index = surf_fit_marker_dof_index[i]; - (*surf_fit_gf)[dof_index] = surf_fit_gf_int(i); - } - - surf_fit_eval_bg_grad->ComputeAtNewPosition( - new_x_sorted, surf_fit_grad_int, ordering); - // Assumes surf_fit_grad and surf_fit_gf share the same space - const int grad_dim = surf_fit_grad->VectorDim(); - const int grad_cnt = surf_fit_grad->Size()/grad_dim; - if (surf_fit_grad->FESpace()->GetOrdering() == Ordering::byNODES) - { - for (int d = 0; d < grad_dim; d++) - { - for (int i = 0; i < cnt; i++) - { - int dof_index = surf_fit_marker_dof_index[i]; - (*surf_fit_grad)[dof_index + d*grad_cnt] = - surf_fit_grad_int(i + d*cnt); - } - } - } - else - { - for (int i = 0; i < cnt; i++) - { - int dof_index = surf_fit_marker_dof_index[i]; - for (int d = 0; d < grad_dim; d++) - { - (*surf_fit_grad)[dof_index*dim + d] = - surf_fit_grad_int(i*dim + d); - } - } - } - - surf_fit_eval_bg_hess->ComputeAtNewPosition( - new_x_sorted, surf_fit_hess_int, ordering); - // Assumes surf_fit_hess and surf_fit_gf share the same space - const int hess_dim = surf_fit_hess->VectorDim(); - const int hess_cnt = surf_fit_hess->Size()/hess_dim; - if (surf_fit_hess->FESpace()->GetOrdering() == Ordering::byNODES) - { - for (int d = 0; d < hess_dim; d++) - { - for (int i = 0; i < cnt; i++) - { - int dof_index = surf_fit_marker_dof_index[i]; - (*surf_fit_hess)[dof_index + d*hess_cnt] = - surf_fit_hess_int(i + d*cnt); - } - } - } - else - { - for (int i = 0; i < cnt; i++) - { - int dof_index = surf_fit_marker_dof_index[i]; - for (int d = 0; d < hess_dim; d++) - { - (*surf_fit_hess)[dof_index*dim + d] = - surf_fit_hess_int(i*dim + d); - } - } - } - } - else - { - surf_fit_eval->ComputeAtNewPosition(x_new, *surf_fit_gf, ordering); - } + RemapSurfaceFittingLevelSetAtNodes(x_new, ordering); } } diff --git a/fem/tmop.hpp b/fem/tmop.hpp index a8aa0a2710..e5b3d54990 100644 --- a/fem/tmop.hpp +++ b/fem/tmop.hpp @@ -1985,6 +1985,10 @@ protected: real_t ComputeUntanglerMaxMuBarrier(const Vector &x, const FiniteElementSpace &fes); + // Remaps the internal surface fitting gridfunction object at provided + // locations. + void RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, + int new_x_ordering); public: /** @param[in] m TMOP_QualityMetric for r-adaptivity (not owned). @param[in] tc Target-matrix construction algorithm to use (not owned). @@ -2105,7 +2109,9 @@ public: /// Parallel support for surface fitting to the zero level set of a function. void EnableSurfaceFitting(const ParGridFunction &s0, const Array &smarker, Coefficient &coeff, - AdaptivityEvaluator &ae); + AdaptivityEvaluator &ae, + AdaptivityEvaluator *aegrad = NULL, + AdaptivityEvaluator *aehess = NULL); /** @brief Fitting of certain DOFs in the current mesh to the zero level set of a function defined on another (finer) source mesh. @@ -2246,6 +2252,9 @@ public: /// Get the surface fitting weight. real_t GetSurfaceFittingWeight(); + /// Copies current surface fitting level-set to provided gridfunction. + void GetSurfaceFittingLevelSet(GridFunction &s0); + /// Computes quantiles needed for UntangleMetrics. Note that in parallel, /// the ParFiniteElementSpace must be passed as argument for consistency /// across MPI ranks. diff --git a/fem/tmop_tools.cpp b/fem/tmop_tools.cpp index afdc99c60e..e06bb919e0 100644 --- a/fem/tmop_tools.cpp +++ b/fem/tmop_tools.cpp @@ -56,7 +56,9 @@ void AdvectorCG::ComputeAtNewPosition(const Vector &new_nodes, new_field_temp(j) = new_field(i + j*ncomp); } } + ComputeAtNewPositionScalar(new_nodes, new_field_temp); + if (fes_ordering == Ordering::byVDIM) { for (int j = 0; j < pnt_cnt; j++) diff --git a/miniapps/meshing/pmesh-fitting.cpp b/miniapps/meshing/pmesh-fitting.cpp index bebb2f75e6..fcbbe41484 100644 --- a/miniapps/meshing/pmesh-fitting.cpp +++ b/miniapps/meshing/pmesh-fitting.cpp @@ -579,6 +579,23 @@ int main (int argc, char *argv[]) } } + // Unify marker across processor boundary + surf_fit_mat_gf.ExchangeFaceNbrData(); + { + + GroupCommunicator &gcomm = surf_fit_mat_gf.ParFESpace()->GroupComm(); + Array gf_array(surf_fit_mat_gf.GetData(), + surf_fit_mat_gf.Size()); + gcomm.Reduce(gf_array, GroupCommunicator::Max); + gcomm.Bcast(gf_array); + } + surf_fit_mat_gf.ExchangeFaceNbrData(); + + for (int i = 0; i < surf_fit_mat_gf.Size(); i++) + { + surf_fit_marker[i] = surf_fit_mat_gf(i) == 1.0; + } + // Set AdaptivityEvaluators for transferring information from initial // mesh to current mesh as it moves during adaptivity. if (adapt_eval == 0) @@ -590,11 +607,8 @@ int main (int argc, char *argv[]) { #ifdef MFEM_USE_GSLIB adapt_surface = new InterpolatorFP; - if (surf_bg_mesh) - { - adapt_grad_surface = new InterpolatorFP; - adapt_hess_surface = new InterpolatorFP; - } + adapt_grad_surface = new InterpolatorFP; + adapt_hess_surface = new InterpolatorFP; #else MFEM_ABORT("MFEM is not built with GSLIB support!"); #endif @@ -604,7 +618,9 @@ int main (int argc, char *argv[]) if (!surf_bg_mesh) { tmop_integ->EnableSurfaceFitting(surf_fit_gf0, surf_fit_marker, - surf_fit_coeff, *adapt_surface); + surf_fit_coeff, *adapt_surface, + adapt_grad_surface, + adapt_hess_surface); } else { @@ -842,9 +858,12 @@ int main (int argc, char *argv[]) if (surface_fit_const > 0.0) { + tmop_integ->GetSurfaceFittingLevelSet(surf_fit_gf0); if (visualization) { - socketstream vis2, vis3; + socketstream vis1, vis2, vis3; + common::VisualizeField(vis1, "localhost", 19916, surf_fit_gf0, + "Level Set", 000, 400, 300, 300); common::VisualizeField(vis2, "localhost", 19916, mat, "Materials", 300, 400, 300, 300); common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf, From 6f3dc3e187a5267e12a002cf625572695ead5214 Mon Sep 17 00:00:00 2001 From: "Mittal, Ketan" Date: Wed, 24 Apr 2024 13:30:35 -0700 Subject: [PATCH 2/9] minor --- fem/tmop_tools.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/fem/tmop_tools.cpp b/fem/tmop_tools.cpp index e06bb919e0..afdc99c60e 100644 --- a/fem/tmop_tools.cpp +++ b/fem/tmop_tools.cpp @@ -56,9 +56,7 @@ void AdvectorCG::ComputeAtNewPosition(const Vector &new_nodes, new_field_temp(j) = new_field(i + j*ncomp); } } - ComputeAtNewPositionScalar(new_nodes, new_field_temp); - if (fes_ordering == Ordering::byVDIM) { for (int j = 0; j < pnt_cnt; j++) From e358c400abdce9dddd9476fb2f8332071ef779b9 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Mon, 20 May 2024 14:48:47 -0700 Subject: [PATCH 3/9] minor --- fem/tmop.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/fem/tmop.cpp b/fem/tmop.cpp index 03df678f49..99d77966dd 100644 --- a/fem/tmop.cpp +++ b/fem/tmop.cpp @@ -3014,41 +3014,41 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0, if (!aegrad) { return; } - // Define spaces for grad and hessian - delete surf_fit_grad; + MFEM_VERIFY(aehess, "AdaptivityEvaluator for Hessians must be provided too."); ParFiniteElementSpace *fes = s0.ParFESpace(); const H1_FECollection *fec = dynamic_cast (surf_fit_gf->FESpace()->FEColl()); if (!fec) { return; } + + // FE space for gradients. + delete surf_fit_grad; H1_FECollection *fec_grad = new H1_FECollection(fec->GetOrder(), dim, fec->GetBasisType()); ParFiniteElementSpace *fes_grad = new ParFiniteElementSpace(pmesh, fec_grad, dim); + // Initial gradients. surf_fit_grad = new GridFunction(fes_grad); surf_fit_grad->MakeOwner(fec_grad); - for (int d = 0; d < dim; d++) { ParGridFunction surf_fit_grad_comp(fes, surf_fit_grad->GetData()+d*s0.Size()); s0.GetDerivative(1, d, surf_fit_grad_comp); } - surf_fit_eval_bg_grad = aegrad; surf_fit_eval_bg_grad->SetParMetaInfo(*pmesh, *fes_grad); surf_fit_eval_bg_grad->SetInitialField(*pmesh->GetNodes(), *surf_fit_grad); - MFEM_VERIFY(aehess,"Specify an adaptivity evaluator for the Hessian terms" - "as well."); + // FE space for Hessians. delete surf_fit_hess; H1_FECollection *fec_hess = new H1_FECollection(fec->GetOrder(), dim, fec->GetBasisType()); ParFiniteElementSpace *fes_hess = new ParFiniteElementSpace(pmesh, fec_hess, dim*dim); + // Initial Hessians. surf_fit_hess = new GridFunction(fes_hess); surf_fit_hess->MakeOwner(fec_hess); - int id = 0; for (int d = 0; d < dim; d++) { @@ -3062,7 +3062,6 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0, id++; } } - surf_fit_eval_bg_hess = aehess; surf_fit_eval_bg_hess->SetParMetaInfo(*pmesh, *fes_hess); surf_fit_eval_bg_hess->SetInitialField(*pmesh->GetNodes(), *surf_fit_hess); @@ -4464,12 +4463,13 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, int new_x_ordering) { if (!surf_fit_gf) { return; } + if (surf_fit_gf_bg && dynamic_cast(surf_fit_eval) && dynamic_cast(surf_fit_eval_bg_grad) && dynamic_cast(surf_fit_eval_bg_hess)) { - // Interpolate information for only DOFs marked for fitting. + // Interpolate information only at DOFs marked for fitting. const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension(); const int cnt = surf_fit_marker_dof_index.Size(); const int total_cnt = new_x.Size()/dim; @@ -4497,7 +4497,7 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, } } - + // Interpolate values of the LS. Vector surf_fit_gf_int, surf_fit_grad_int, surf_fit_hess_int; surf_fit_eval->ComputeAtNewPosition(new_x_sorted, surf_fit_gf_int, new_x_ordering); @@ -4507,10 +4507,9 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, (*surf_fit_gf)[dof_index] = surf_fit_gf_int(i); } - + // Interpolate gradients of the LS. surf_fit_eval_bg_grad->ComputeAtNewPosition(new_x_sorted, surf_fit_grad_int, new_x_ordering); - // Assumes surf_fit_grad and surf_fit_gf share the same space const int grad_dim = surf_fit_grad->VectorDim(); const int grad_cnt = surf_fit_grad->Size()/grad_dim; @@ -4539,6 +4538,7 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, } } + // Interpolate Hessians of the LS. surf_fit_eval_bg_hess->ComputeAtNewPosition(new_x_sorted, surf_fit_hess_int, new_x_ordering); // Assumes surf_fit_hess and surf_fit_gf share the same space From 9cfae52d1ef8987fbd398cddb8ff39ebdcf8817c Mon Sep 17 00:00:00 2001 From: "Mittal, Ketan" Date: Tue, 21 May 2024 10:07:48 -0700 Subject: [PATCH 4/9] minor changes based on reviewer comments --- fem/tmop.cpp | 30 ++++++++++++++---------------- fem/tmop.hpp | 3 --- miniapps/meshing/pmesh-fitting.cpp | 3 ++- 3 files changed, 16 insertions(+), 20 deletions(-) diff --git a/fem/tmop.cpp b/fem/tmop.cpp index 99d77966dd..bf107639a7 100644 --- a/fem/tmop.cpp +++ b/fem/tmop.cpp @@ -2999,6 +2999,11 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0, ParMesh *pmesh = s0.ParFESpace()->GetParMesh(); MFEM_VERIFY(pmesh->GetNodes()->Size() == dim*s0.Size(), "Mesh and level-set polynomial order must be the same."); + const H1_FECollection *fec = dynamic_cast + (s0.FESpace()->FEColl()); + MFEM_VERIFY(fec, "Only H1_FECollection is supported for the surface fitting " + "grid function."); + delete surf_fit_gf; surf_fit_gf = new GridFunction(s0); @@ -3018,10 +3023,6 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0, ParFiniteElementSpace *fes = s0.ParFESpace(); - const H1_FECollection *fec = dynamic_cast - (surf_fit_gf->FESpace()->FEColl()); - if (!fec) { return; } - // FE space for gradients. delete surf_fit_grad; H1_FECollection *fec_grad = new H1_FECollection(fec->GetOrder(), dim, @@ -4319,15 +4320,6 @@ real_t TMOP_Integrator::GetSurfaceFittingWeight() return 0.0; } -void TMOP_Integrator::GetSurfaceFittingLevelSet(GridFunction &s0) -{ - - MFEM_VERIFY(s0.Size() == surf_fit_gf->Size(), - "Provided function and internal level-set function should" - "be of same size."); - s0 = *surf_fit_gf; -} - void TMOP_Integrator::EnableNormalization(const GridFunction &x) { ComputeNormalizationEnergies(x, metric_normal, lim_normal, surf_fit_normal); @@ -4463,11 +4455,17 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, int new_x_ordering) { if (!surf_fit_gf) { return; } - - if (surf_fit_gf_bg && - dynamic_cast(surf_fit_eval) && + bool optimized_remap = false; // only remap at nodes for fitting. +#ifdef MFEM_USE_GSLIB + if (dynamic_cast(surf_fit_eval) && dynamic_cast(surf_fit_eval_bg_grad) && dynamic_cast(surf_fit_eval_bg_hess)) + { + optimized_remap = true; + } +#endif + + if (surf_fit_gf_bg && optimized_remap) { // Interpolate information only at DOFs marked for fitting. const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension(); diff --git a/fem/tmop.hpp b/fem/tmop.hpp index e5b3d54990..3b940dc51e 100644 --- a/fem/tmop.hpp +++ b/fem/tmop.hpp @@ -2252,9 +2252,6 @@ public: /// Get the surface fitting weight. real_t GetSurfaceFittingWeight(); - /// Copies current surface fitting level-set to provided gridfunction. - void GetSurfaceFittingLevelSet(GridFunction &s0); - /// Computes quantiles needed for UntangleMetrics. Note that in parallel, /// the ParFiniteElementSpace must be passed as argument for consistency /// across MPI ranks. diff --git a/miniapps/meshing/pmesh-fitting.cpp b/miniapps/meshing/pmesh-fitting.cpp index fcbbe41484..34ae27f004 100644 --- a/miniapps/meshing/pmesh-fitting.cpp +++ b/miniapps/meshing/pmesh-fitting.cpp @@ -858,7 +858,8 @@ int main (int argc, char *argv[]) if (surface_fit_const > 0.0) { - tmop_integ->GetSurfaceFittingLevelSet(surf_fit_gf0); + adapt_surface->ComputeAtNewPosition(x, surf_fit_gf0, + x.FESpace()->GetOrdering()); if (visualization) { socketstream vis1, vis2, vis3; From a910f49710f11ec28baf04575baf839e89c487e5 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Tue, 21 May 2024 14:01:57 -0700 Subject: [PATCH 5/9] minor --- fem/tmop.cpp | 7 ++++--- miniapps/meshing/pmesh-fitting.cpp | 1 - 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/fem/tmop.cpp b/fem/tmop.cpp index bf107639a7..bfefd8240b 100644 --- a/fem/tmop.cpp +++ b/fem/tmop.cpp @@ -4455,17 +4455,18 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, int new_x_ordering) { if (!surf_fit_gf) { return; } - bool optimized_remap = false; // only remap at nodes for fitting. + + bool remap_only_fitted_nodes = false; // only remap at nodes for fitting. #ifdef MFEM_USE_GSLIB if (dynamic_cast(surf_fit_eval) && dynamic_cast(surf_fit_eval_bg_grad) && dynamic_cast(surf_fit_eval_bg_hess)) { - optimized_remap = true; + remap_only_fitted_nodes = true; } #endif - if (surf_fit_gf_bg && optimized_remap) + if (surf_fit_gf_bg && remap_only_fitted_nodes) { // Interpolate information only at DOFs marked for fitting. const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension(); diff --git a/miniapps/meshing/pmesh-fitting.cpp b/miniapps/meshing/pmesh-fitting.cpp index 34ae27f004..b40d448f21 100644 --- a/miniapps/meshing/pmesh-fitting.cpp +++ b/miniapps/meshing/pmesh-fitting.cpp @@ -582,7 +582,6 @@ int main (int argc, char *argv[]) // Unify marker across processor boundary surf_fit_mat_gf.ExchangeFaceNbrData(); { - GroupCommunicator &gcomm = surf_fit_mat_gf.ParFESpace()->GroupComm(); Array gf_array(surf_fit_mat_gf.GetData(), surf_fit_mat_gf.Size()); From e195a709ffe40f41dd6d2d1d0e7b3b5b236a85d2 Mon Sep 17 00:00:00 2001 From: "Mittal, Ketan" Date: Tue, 21 May 2024 15:31:22 -0700 Subject: [PATCH 6/9] better documentation and changed some logic --- fem/tmop.cpp | 99 +++++++++++++++++++++++++++------------------------- fem/tmop.hpp | 18 ++++++---- 2 files changed, 63 insertions(+), 54 deletions(-) diff --git a/fem/tmop.cpp b/fem/tmop.cpp index bfefd8240b..f78ff20ef3 100644 --- a/fem/tmop.cpp +++ b/fem/tmop.cpp @@ -2949,6 +2949,15 @@ void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &s0, MFEM_VERIFY(surf_fit_pos == NULL, "Using both fitting approaches is not supported."); + const int dim = s0.FESpace()->GetMesh()->Dimension(); + Mesh *mesh = s0.FESpace()->GetMesh(); + MFEM_VERIFY(mesh->GetNodes()->Size() == dim*s0.Size(), + "Mesh and level-set polynomial order must be the same."); + const H1_FECollection *fec = dynamic_cast + (s0.FESpace()->FEColl()); + MFEM_VERIFY(fec, "Only H1_FECollection is supported for the surface fitting " + "grid function."); + delete surf_fit_gf; surf_fit_gf = new GridFunction(s0); surf_fit_gf->CountElementsPerVDof(surf_fit_dof_count); @@ -3015,7 +3024,6 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0, surf_fit_eval->SetParMetaInfo(*pmesh, *s0.ParFESpace()); surf_fit_eval->SetInitialField (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf); - surf_fit_gf_bg = false; if (!aegrad) { return; } @@ -3037,9 +3045,9 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0, ParGridFunction surf_fit_grad_comp(fes, surf_fit_grad->GetData()+d*s0.Size()); s0.GetDerivative(1, d, surf_fit_grad_comp); } - surf_fit_eval_bg_grad = aegrad; - surf_fit_eval_bg_grad->SetParMetaInfo(*pmesh, *fes_grad); - surf_fit_eval_bg_grad->SetInitialField(*pmesh->GetNodes(), *surf_fit_grad); + surf_fit_eval_grad = aegrad; + surf_fit_eval_grad->SetParMetaInfo(*pmesh, *fes_grad); + surf_fit_eval_grad->SetInitialField(*pmesh->GetNodes(), *surf_fit_grad); // FE space for Hessians. delete surf_fit_hess; @@ -3063,21 +3071,27 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0, id++; } } - surf_fit_eval_bg_hess = aehess; - surf_fit_eval_bg_hess->SetParMetaInfo(*pmesh, *fes_hess); - surf_fit_eval_bg_hess->SetInitialField(*pmesh->GetNodes(), *surf_fit_hess); + surf_fit_eval_hess = aehess; + surf_fit_eval_hess->SetParMetaInfo(*pmesh, *fes_hess); + surf_fit_eval_hess->SetInitialField(*pmesh->GetNodes(), *surf_fit_hess); - surf_fit_gf_bg = true; // Store DOF indices that are marked for fitting. Used to reduce work for // transferring information between source/background and current mesh. surf_fit_marker_dof_index.SetSize(0); - for (int i = 0; i < surf_fit_marker->Size(); i++) +#ifdef MFEM_USE_GSLIB + if (dynamic_cast(surf_fit_eval) && + dynamic_cast(surf_fit_eval_grad) && + dynamic_cast(surf_fit_eval_hess)) { - if ((*surf_fit_marker)[i] == true) + for (int i = 0; i < surf_fit_marker->Size(); i++) { - surf_fit_marker_dof_index.Append(i); + if ((*surf_fit_marker)[i] == true) + { + surf_fit_marker_dof_index.Append(i); + } } } +#endif *surf_fit_grad = 0.0; *surf_fit_hess = 0.0; @@ -3116,11 +3130,11 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource( delete surf_fit_grad; surf_fit_grad = new GridFunction(s0_grad); *surf_fit_grad = 0.0; - surf_fit_eval_bg_grad = &age; - surf_fit_eval_bg_hess = &ahe; - surf_fit_eval_bg_grad->SetParMetaInfo(*s_bg_grad.ParFESpace()->GetParMesh(), - *s_bg_grad.ParFESpace()); - surf_fit_eval_bg_grad->SetInitialField + surf_fit_eval_grad = &age; + surf_fit_eval_hess = &ahe; + surf_fit_eval_grad->SetParMetaInfo(*s_bg_grad.ParFESpace()->GetParMesh(), + *s_bg_grad.ParFESpace()); + surf_fit_eval_grad->SetInitialField (*s_bg_grad.FESpace()->GetMesh()->GetNodes(), s_bg_grad); // Setup for Hessian on background mesh @@ -3131,9 +3145,9 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource( delete surf_fit_hess; surf_fit_hess = new GridFunction(s0_hess); *surf_fit_hess = 0.0; - surf_fit_eval_bg_hess->SetParMetaInfo(*s_bg_hess.ParFESpace()->GetParMesh(), - *s_bg_hess.ParFESpace()); - surf_fit_eval_bg_hess->SetInitialField + surf_fit_eval_hess->SetParMetaInfo(*s_bg_hess.ParFESpace()->GetParMesh(), + *s_bg_hess.ParFESpace()); + surf_fit_eval_hess->SetInitialField (*s_bg_hess.FESpace()->GetMesh()->GetNodes(), s_bg_hess); // Count number of zones that share each of the DOFs @@ -3938,7 +3952,7 @@ void TMOP_Integrator::AssembleElemVecSurfFit(const FiniteElement &el_x, Vector sigma_e(dof_s); DenseMatrix surf_fit_grad_e(dof_s, dim); - if (surf_fit_gf || surf_fit_gf_bg) + if (surf_fit_gf) { surf_fit_gf->GetSubVector(vdofs, sigma_e); @@ -3946,7 +3960,7 @@ void TMOP_Integrator::AssembleElemVecSurfFit(const FiniteElement &el_x, // The FE coefficients of the gradient go in surf_fit_grad_e. Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim); DenseMatrix grad_phys; // This will be (dof x dim, dof). - if (surf_fit_gf_bg) + if (surf_fit_grad) { surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs); surf_fit_grad->GetSubVector(dofs, grad_ptr); @@ -4020,7 +4034,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x, Vector sigma_e(dof_s); DenseMatrix surf_fit_grad_e(dof_s, dim); DenseMatrix surf_fit_hess_e(dof_s, dim*dim); - if (surf_fit_gf || surf_fit_gf_bg) + if (surf_fit_gf) { surf_fit_gf->GetSubVector(vdofs, sigma_e); @@ -4028,7 +4042,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x, // The FE coefficients of the gradient go in surf_fit_grad_e. Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim); DenseMatrix grad_phys; // This will be (dof x dim, dof). - if (surf_fit_gf_bg) + if (surf_fit_grad) { surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs); surf_fit_grad->GetSubVector(dofs, grad_ptr); @@ -4042,7 +4056,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x, // Project the Hessian of sigma in the same space. // The FE coefficients of the Hessian go in surf_fit_hess_e. Vector hess_ptr(surf_fit_hess_e.GetData(), dof_s*dim*dim); - if (surf_fit_gf_bg) + if (surf_fit_hess) { surf_fit_hess->FESpace()->GetElementVDofs(el_id, dofs); surf_fit_hess->GetSubVector(dofs, hess_ptr); @@ -4069,7 +4083,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x, Tpr.SetIntPoint(&ip); real_t w = surf_fit_normal * surf_fit_coeff->Eval(Tpr, ip); - if (surf_fit_gf || surf_fit_gf_bg) + if (surf_fit_gf) { Vector gg_ptr(surf_fit_hess_s.GetData(), dim * dim); surf_fit_hess_e.GetRow(s, gg_ptr); @@ -4456,17 +4470,7 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, { if (!surf_fit_gf) { return; } - bool remap_only_fitted_nodes = false; // only remap at nodes for fitting. -#ifdef MFEM_USE_GSLIB - if (dynamic_cast(surf_fit_eval) && - dynamic_cast(surf_fit_eval_bg_grad) && - dynamic_cast(surf_fit_eval_bg_hess)) - { - remap_only_fitted_nodes = true; - } -#endif - - if (surf_fit_gf_bg && remap_only_fitted_nodes) + if (surf_fit_marker_dof_index.Size()) { // Interpolate information only at DOFs marked for fitting. const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension(); @@ -4507,8 +4511,8 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, } // Interpolate gradients of the LS. - surf_fit_eval_bg_grad->ComputeAtNewPosition(new_x_sorted, surf_fit_grad_int, - new_x_ordering); + surf_fit_eval_grad->ComputeAtNewPosition(new_x_sorted, surf_fit_grad_int, + new_x_ordering); // Assumes surf_fit_grad and surf_fit_gf share the same space const int grad_dim = surf_fit_grad->VectorDim(); const int grad_cnt = surf_fit_grad->Size()/grad_dim; @@ -4538,8 +4542,8 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, } // Interpolate Hessians of the LS. - surf_fit_eval_bg_hess->ComputeAtNewPosition(new_x_sorted, surf_fit_hess_int, - new_x_ordering); + surf_fit_eval_hess->ComputeAtNewPosition(new_x_sorted, surf_fit_hess_int, + new_x_ordering); // Assumes surf_fit_hess and surf_fit_gf share the same space const int hess_dim = surf_fit_hess->VectorDim(); const int hess_cnt = surf_fit_hess->Size()/hess_dim; @@ -4572,15 +4576,15 @@ void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x, else { surf_fit_eval->ComputeAtNewPosition(new_x, *surf_fit_gf, new_x_ordering); - if (surf_fit_eval_bg_grad) + if (surf_fit_eval_grad) { - surf_fit_eval_bg_grad->ComputeAtNewPosition(new_x, *surf_fit_grad, - new_x_ordering); + surf_fit_eval_grad->ComputeAtNewPosition(new_x, *surf_fit_grad, + new_x_ordering); } - if (surf_fit_eval_bg_hess) + if (surf_fit_eval_hess) { - surf_fit_eval_bg_hess->ComputeAtNewPosition(new_x, *surf_fit_hess, - new_x_ordering); + surf_fit_eval_hess->ComputeAtNewPosition(new_x, *surf_fit_hess, + new_x_ordering); } } } @@ -4615,7 +4619,8 @@ UpdateAfterMeshPositionChange(const Vector &x_new, adapt_lim_eval->ComputeAtNewPosition(x_new, *adapt_lim_gf, ordering); } - // Update surf_fit_gf if surface fitting is enabled. + // Update surf_fit_gf (and optionally its gradients) if surface + // fitting is enabled. if (surf_fit_gf) { RemapSurfaceFittingLevelSetAtNodes(x_new, ordering); diff --git a/fem/tmop.hpp b/fem/tmop.hpp index 3b940dc51e..9d929bcaa7 100644 --- a/fem/tmop.hpp +++ b/fem/tmop.hpp @@ -1784,12 +1784,12 @@ protected: // Fitting to given physical positions. TMOP_QuadraticLimiter *surf_fit_limiter; // Owned. Created internally. const GridFunction *surf_fit_pos; // Not owned. Positions to fit. - real_t surf_fit_normal; - bool surf_fit_gf_bg; - GridFunction *surf_fit_grad, *surf_fit_hess; - AdaptivityEvaluator *surf_fit_eval_bg_grad, *surf_fit_eval_bg_hess; - Array surf_fit_dof_count; - Array surf_fit_marker_dof_index; + real_t surf_fit_normal; // Normalization factor. + bool surf_fit_gf_bg; // Use background mesh for fitting. + GridFunction *surf_fit_grad, *surf_fit_hess; // Owned. Created internally. + AdaptivityEvaluator *surf_fit_eval_grad, *surf_fit_eval_hess; // Not owned. + Array surf_fit_dof_count; // Number of dofs per node. + Array surf_fit_marker_dof_index; // Indices of nodes to fit. DiscreteAdaptTC *discr_tc; @@ -2006,7 +2006,7 @@ public: surf_fit_limiter(NULL), surf_fit_pos(NULL), surf_fit_normal(1.0), surf_fit_gf_bg(false), surf_fit_grad(NULL), surf_fit_hess(NULL), - surf_fit_eval_bg_grad(NULL), surf_fit_eval_bg_hess(NULL), + surf_fit_eval_grad(NULL), surf_fit_eval_hess(NULL), discr_tc(dynamic_cast(tc)), fdflag(false), dxscale(1.0e3), fd_call_flag(false), exact_action(false) { PA.enabled = false; } @@ -2107,6 +2107,10 @@ public: #ifdef MFEM_USE_MPI /// Parallel support for surface fitting to the zero level set of a function. + /// Here, we add two optional inputs: @a aegrad and @a aehess. When provided, + /// the first and second derivative of the input level set are computed on + /// the initial mesh, and @a aegrad and @a aehess are used to remap grad_s(x) + /// from grad_s0(x0) and hess_s(x) from hess_s0(x0), respectively. void EnableSurfaceFitting(const ParGridFunction &s0, const Array &smarker, Coefficient &coeff, AdaptivityEvaluator &ae, From 69a4a38053d0323cd656facfdbc8d4e724148220 Mon Sep 17 00:00:00 2001 From: "Mittal, Ketan" Date: Tue, 21 May 2024 15:36:49 -0700 Subject: [PATCH 7/9] remove unneeded flag --- fem/tmop.cpp | 1 - fem/tmop.hpp | 4 +--- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/fem/tmop.cpp b/fem/tmop.cpp index f78ff20ef3..78760c6af3 100644 --- a/fem/tmop.cpp +++ b/fem/tmop.cpp @@ -3116,7 +3116,6 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource( surf_fit_coeff = &coeff; surf_fit_eval = &ae; - surf_fit_gf_bg = true; surf_fit_eval->SetParMetaInfo(*s_bg.ParFESpace()->GetParMesh(), *s_bg.ParFESpace()); surf_fit_eval->SetInitialField diff --git a/fem/tmop.hpp b/fem/tmop.hpp index 9d929bcaa7..841637f1ee 100644 --- a/fem/tmop.hpp +++ b/fem/tmop.hpp @@ -1785,7 +1785,6 @@ protected: TMOP_QuadraticLimiter *surf_fit_limiter; // Owned. Created internally. const GridFunction *surf_fit_pos; // Not owned. Positions to fit. real_t surf_fit_normal; // Normalization factor. - bool surf_fit_gf_bg; // Use background mesh for fitting. GridFunction *surf_fit_grad, *surf_fit_hess; // Owned. Created internally. AdaptivityEvaluator *surf_fit_eval_grad, *surf_fit_eval_hess; // Not owned. Array surf_fit_dof_count; // Number of dofs per node. @@ -2004,8 +2003,7 @@ public: surf_fit_marker(NULL), surf_fit_coeff(NULL), surf_fit_gf(NULL), surf_fit_eval(NULL), surf_fit_limiter(NULL), surf_fit_pos(NULL), - surf_fit_normal(1.0), - surf_fit_gf_bg(false), surf_fit_grad(NULL), surf_fit_hess(NULL), + surf_fit_normal(1.0), surf_fit_grad(NULL), surf_fit_hess(NULL), surf_fit_eval_grad(NULL), surf_fit_eval_hess(NULL), discr_tc(dynamic_cast(tc)), fdflag(false), dxscale(1.0e3), fd_call_flag(false), exact_action(false) From 4ee1bcd56181e3bdf1b3fed70f1609adb5ebb05b Mon Sep 17 00:00:00 2001 From: "Mittal, Ketan" Date: Tue, 21 May 2024 16:02:38 -0700 Subject: [PATCH 8/9] double -> real_t --- miniapps/meshing/pmesh-fitting.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/miniapps/meshing/pmesh-fitting.cpp b/miniapps/meshing/pmesh-fitting.cpp index b40d448f21..6395ba92e7 100644 --- a/miniapps/meshing/pmesh-fitting.cpp +++ b/miniapps/meshing/pmesh-fitting.cpp @@ -583,9 +583,9 @@ int main (int argc, char *argv[]) surf_fit_mat_gf.ExchangeFaceNbrData(); { GroupCommunicator &gcomm = surf_fit_mat_gf.ParFESpace()->GroupComm(); - Array gf_array(surf_fit_mat_gf.GetData(), + Array gf_array(surf_fit_mat_gf.GetData(), surf_fit_mat_gf.Size()); - gcomm.Reduce(gf_array, GroupCommunicator::Max); + gcomm.Reduce(gf_array, GroupCommunicator::Max); gcomm.Bcast(gf_array); } surf_fit_mat_gf.ExchangeFaceNbrData(); From 1f9c75585ee4f6957dc77df11b9e6f2db3e6b161 Mon Sep 17 00:00:00 2001 From: "Mittal, Ketan" Date: Wed, 22 May 2024 11:59:42 -0700 Subject: [PATCH 9/9] add mfem_use_mpi guard in gslib --- fem/gslib.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/fem/gslib.cpp b/fem/gslib.cpp index eb040c25e7..7b70e27f38 100644 --- a/fem/gslib.cpp +++ b/fem/gslib.cpp @@ -1393,8 +1393,10 @@ GSOPGSLIB::~GSOPGSLIB() void GSOPGSLIB::UpdateIdentifiers(const Array &ids) { long long minval = ids.Min(); +#ifdef MFEM_USE_MPI MPI_Allreduce(MPI_IN_PLACE, &minval, 1, MPI_LONG_LONG_INT, MPI_MIN, gsl_comm->c); +#endif MFEM_VERIFY(minval >= 0, "Unique identifier cannot be negative."); if (gsl_data != NULL) { gslib_gs_free(gsl_data); } num_ids = ids.Size();