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/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(); diff --git a/fem/pgridfunc.cpp b/fem/pgridfunc.cpp index 3d686d34dc..282ec9aa1d 100644 --- a/fem/pgridfunc.cpp +++ b/fem/pgridfunc.cpp @@ -520,7 +520,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 cbb353c073..d5e388132b 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 719b261f3d..92277376db 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); @@ -2987,12 +2996,24 @@ 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."); + 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); s0.CountElementsPerVDof(surf_fit_dof_count); @@ -3000,11 +3021,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; } + + MFEM_VERIFY(aehess, "AdaptivityEvaluator for Hessians must be provided too."); + + ParFiniteElementSpace *fes = s0.ParFESpace(); + + // 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_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; + 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++) + { + 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_hess = aehess; + surf_fit_eval_hess->SetParMetaInfo(*pmesh, *fes_hess); + surf_fit_eval_hess->SetInitialField(*pmesh->GetNodes(), *surf_fit_hess); + + // 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); +#ifdef MFEM_USE_GSLIB + if (dynamic_cast(surf_fit_eval) && + dynamic_cast(surf_fit_eval_grad) && + dynamic_cast(surf_fit_eval_hess)) + { + for (int i = 0; i < surf_fit_marker->Size(); 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; } void TMOP_Integrator::EnableSurfaceFittingFromSource( @@ -3022,12 +3112,10 @@ 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; - surf_fit_gf_bg = true; surf_fit_eval->SetParMetaInfo(*s_bg.ParFESpace()->GetParMesh(), *s_bg.ParFESpace()); surf_fit_eval->SetInitialField @@ -3044,11 +3132,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 @@ -3059,9 +3147,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 @@ -3866,7 +3954,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); @@ -3874,7 +3962,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); @@ -3948,7 +4036,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); @@ -3956,7 +4044,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); @@ -3970,7 +4058,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); @@ -3997,7 +4085,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); @@ -4379,6 +4467,130 @@ 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_marker_dof_index.Size()) + { + // 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; + 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); + } + } + } + + // 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); + 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); + } + + // Interpolate gradients of the LS. + 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; + 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); + } + } + } + + // Interpolate Hessians of the LS. + 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; + 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_grad) + { + surf_fit_eval_grad->ComputeAtNewPosition(new_x, *surf_fit_grad, + new_x_ordering); + } + if (surf_fit_eval_hess) + { + surf_fit_eval_hess->ComputeAtNewPosition(new_x, *surf_fit_hess, + new_x_ordering); + } + } +} + void TMOP_Integrator:: UpdateAfterMeshPositionChange(const Vector &x_new, const FiniteElementSpace &x_fes) @@ -4409,112 +4621,11 @@ 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) { - 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..841637f1ee 100644 --- a/fem/tmop.hpp +++ b/fem/tmop.hpp @@ -1784,12 +1784,11 @@ 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. + 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; @@ -1985,6 +1984,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). @@ -2000,9 +2003,8 @@ 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_eval_bg_grad(NULL), surf_fit_eval_bg_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) { PA.enabled = false; } @@ -2103,9 +2105,15 @@ 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); + 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. diff --git a/miniapps/meshing/pmesh-fitting.cpp b/miniapps/meshing/pmesh-fitting.cpp index 93a00e3084..08db3f74a6 100644 --- a/miniapps/meshing/pmesh-fitting.cpp +++ b/miniapps/meshing/pmesh-fitting.cpp @@ -580,6 +580,22 @@ 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) @@ -591,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 @@ -605,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 { @@ -837,9 +852,13 @@ int main (int argc, char *argv[]) if (surface_fit_const > 0.0) { + adapt_surface->ComputeAtNewPosition(x, surf_fit_gf0, + x.FESpace()->GetOrdering()); 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,