diff --git a/cases/condif2/condif.c b/cases/condif2/condif.c new file mode 100644 index 0000000..3d9b4a3 --- /dev/null +++ b/cases/condif2/condif.c @@ -0,0 +1,35 @@ +//-------------------------------------------------------------- +// Solution function for convection-diffusion +// +// To compile run, for example: +// mpicc -shared -o libfun.so -fPIC condif.c +// gcc -shared -o libfun.so -fPIC condif.c +//-------------------------------------------------------------- + +#include +#include + + +double sol_phi(double *coord, int dim, double time) +{ + return 0.0; +} + +void advection(double *coord, int dim, double time, double *adv, int vdim) +{ + double r = sqrt(coord[0]*coord[0] + coord[1]*coord[1]); + + adv[0] = -coord[0]/r;//-sqrt(2.0)/2.0; + adv[1] = -coord[1]/r;//-sqrt(2.0)/2.0; +} + +double force(double *coord, int dim, double time) +{ + return 1.0; +} + +double mu(double *coord, int dim, double time) +{ + return 0.00001; +} + diff --git a/cases/condif2/condif_ms.c b/cases/condif2/condif_ms.c new file mode 100644 index 0000000..8188066 --- /dev/null +++ b/cases/condif2/condif_ms.c @@ -0,0 +1,56 @@ +//-------------------------------------------------------------- +// Solution function for convection-diffusion -- manufactured solution +// +// To compile run, for example: +// mpicc -shared -o libfun.so -fPIC condif_ms.c +// gcc -shared -o libfun.so -fPIC condif_ms.c +//-------------------------------------------------------------- + +# TO BE MODIFIED ??? +bla + +#include +#include + +void advection(double *x, int dim, double t, double *adv, int vdim) +{ + adv[0] = sqrt(2.0)/2.0; + adv[1] = sqrt(2.0)/2.0; +} + +double mu(double *x, int dim, double t) +{ + return 0.01; +} + +double sol_phi(double *x, int dim, double t) +{ + return sin(M_PI*x[0])*sin(M_PI*x[1]); +} + +void grad_phi(double *x, int dim, double t, double *grad, int vdim) +{ + grad[0] = M_PI*cos(M_PI*x[0])*sin(M_PI*x[1]); + grad[1] = M_PI*sin(M_PI*x[0])*cos(M_PI*x[1]); +} + +double lap_phi(double *x, int dim, double t) +{ + return -2*M_PI*M_PI*sin(M_PI*x[0])*sin(M_PI*x[1]); +} + +double force(double *x, int dim, double t) +{ + double a[dim]; + double grad[dim]; + advection(x, dim, t, a, dim); + grad_phi(x, dim, t, grad, dim); + + double res = - mu(x, dim, t)*lap_phi(x, dim, t); + for (int i = 0; i < dim; i++) + { + res += a[i]*grad[i]; + } + + return res; +} diff --git a/cases/condif2/quart-circle-eps.mesh b/cases/condif2/quart-circle-eps.mesh new file mode 100644 index 0000000..f1d311f --- /dev/null +++ b/cases/condif2/quart-circle-eps.mesh @@ -0,0 +1,64 @@ +MFEM NURBS mesh v1.0 + +# +# MFEM Geometry Types (see fem/geom.hpp): +# +# SEGMENT = 1 +# SQUARE = 3 +# CUBE = 5 +# + +dimension +2 + +elements +1 +1 3 0 1 2 3 + +boundary +4 +1 1 0 1 +3 1 2 3 +4 1 3 0 +2 1 1 2 + +edges +4 +0 0 1 +0 3 2 +1 0 3 +1 1 2 + +vertices +4 + +knotvectors +2 +2 3 0 0 0 1 1 1 +2 3 0 0 0 1 1 1 + +weights +1 +1 +1 +1 +0.707106781 +0.707106781 +1 +1 +0.707106781 + +FiniteElementSpace +FiniteElementCollection: NURBS2 +VDim: 2 +Ordering: 1 + +0 0.001 +0.001 0 +1 0 +0 1 +0.001 0.001 +1 1 +0 0.001 +0.001 0 +0.001 0.001 diff --git a/cases/condif2/quart-circle.mesh b/cases/condif2/quart-circle.mesh new file mode 100644 index 0000000..b94ad05 --- /dev/null +++ b/cases/condif2/quart-circle.mesh @@ -0,0 +1,64 @@ +MFEM NURBS mesh v1.0 + +# +# MFEM Geometry Types (see fem/geom.hpp): +# +# SEGMENT = 1 +# SQUARE = 3 +# CUBE = 5 +# + +dimension +2 + +elements +1 +1 3 0 1 2 3 + +boundary +4 +1 1 0 1 +3 1 2 3 +4 1 3 0 +2 1 1 2 + +edges +4 +0 0 1 +0 3 2 +1 0 3 +1 1 2 + +vertices +4 + +knotvectors +2 +2 3 0 0 0 1 1 1 +2 3 0 0 0 1 1 1 + +weights +1 +1 +1 +1 +0.707106781 +0.707106781 +1 +1 +0.707106781 + +FiniteElementSpace +FiniteElementCollection: NURBS2 +VDim: 2 +Ordering: 1 + +0 0.00 +0.00 0 +1 0 +0 1 +0.00 0.00 +1 1 +0 0.5 +0.5 0 +0.5 0.5 diff --git a/cases/condif2/run-nurbs b/cases/condif2/run-nurbs new file mode 100644 index 0000000..b01b745 --- /dev/null +++ b/cases/condif2/run-nurbs @@ -0,0 +1,9 @@ +#!/bin/bash +gcc -shared -o libfun.so -fPIC condif.c + +exe=../../build/supg/supg +mesh=quart-circle.mesh + +rm output* log* + +mpirun -n 1 $exe -m $mesh -o 2 -r 3 --strong-bdr "1 2 3 4" -l libfun.so | tee log0 diff --git a/cases/condif3/condif.c b/cases/condif3/condif.c new file mode 100644 index 0000000..d405d02 --- /dev/null +++ b/cases/condif3/condif.c @@ -0,0 +1,35 @@ +//-------------------------------------------------------------- +// Solution function for convection-diffusion +// +// To compile run, for example: +// mpicc -shared -o libfun.so -fPIC condif.c +// gcc -shared -o libfun.so -fPIC condif.c +//-------------------------------------------------------------- + +#include +#include + + +double sol_phi(double *coord, int dim, double time) +{ + double r = sqrt(coord[0]*coord[0] + coord[1]*coord[1]); + + return 0.5*(1.0 - cos(2*M_PI*fmin(1.0, r))); +} + +void advection(double *coord, int dim, double time, double *adv, int vdim) +{ + adv[0] = coord[1]; + adv[1] = -coord[0]; +} + +double force(double *coord, int dim, double time) +{ + return 0.0; +} + +double mu(double *coord, int dim, double time) +{ + return 0.0; +} + diff --git a/cases/condif3/condif_ms.c b/cases/condif3/condif_ms.c new file mode 100644 index 0000000..57f3080 --- /dev/null +++ b/cases/condif3/condif_ms.c @@ -0,0 +1,55 @@ +//-------------------------------------------------------------- +// Solution function for convection-diffusion -- manufactured solution +// +// To compile run, for example: +// mpicc -shared -o libfun.so -fPIC condif_ms.c +// gcc -shared -o libfun.so -fPIC condif_ms.c +//-------------------------------------------------------------- + +## ?? + +#include +#include + +void advection(double *x, int dim, double t, double *adv, int vdim) +{ + adv[0] = sqrt(2.0)/2.0; + adv[1] = sqrt(2.0)/2.0; +} + +double mu(double *x, int dim, double t) +{ + return 0.01; +} + +double sol_phi(double *x, int dim, double t) +{ + return sin(M_PI*x[0])*sin(M_PI*x[1]); +} + +void grad_phi(double *x, int dim, double t, double *grad, int vdim) +{ + grad[0] = M_PI*cos(M_PI*x[0])*sin(M_PI*x[1]); + grad[1] = M_PI*sin(M_PI*x[0])*cos(M_PI*x[1]); +} + +double lap_phi(double *x, int dim, double t) +{ + return -2*M_PI*M_PI*sin(M_PI*x[0])*sin(M_PI*x[1]); +} + +double force(double *x, int dim, double t) +{ + double a[dim]; + double grad[dim]; + advection(x, dim, t, a, dim); + grad_phi(x, dim, t, grad, dim); + + double res = - mu(x, dim, t)*lap_phi(x, dim, t); + for (int i = 0; i < dim; i++) + { + res += a[i]*grad[i]; + } + + return res; +} diff --git a/cases/condif3/run-nurbs b/cases/condif3/run-nurbs new file mode 100644 index 0000000..6425897 --- /dev/null +++ b/cases/condif3/run-nurbs @@ -0,0 +1,10 @@ +#!/bin/bash +gcc -shared -o libfun.so -fPIC condif.c + +exe=../../build/supg/supg +mesh=square.mesh +mesh=square-slit.mesh + +rm output* log* + +mpirun -n 1 $exe -m $mesh -o 2 -r 3 --strong-bdr "1 2" -l libfun.so | tee log0 diff --git a/cases/condif3/square-slit.mesh b/cases/condif3/square-slit.mesh new file mode 100644 index 0000000..27add35 --- /dev/null +++ b/cases/condif3/square-slit.mesh @@ -0,0 +1,110 @@ +MFEM NURBS mesh v1.0 + +# +# MFEM Geometry Types (see fem/geom.hpp): +# +# SEGMENT = 1 +# SQUARE = 3 +# CUBE = 5 +# + +dimension +2 + +elements +4 +1 3 0 1 5 4 +1 3 2 3 6 5 +1 3 4 5 8 7 +1 3 5 6 9 8 + +boundary +10 +1 1 0 1 +1 1 1 5 +1 1 4 0 +1 1 2 3 +1 1 3 6 +2 1 5 2 +1 1 8 7 +1 1 7 4 +1 1 6 9 +1 1 9 8 + +edges +13 +0 0 1 +1 1 5 +0 4 5 +1 0 4 +2 2 3 +3 3 6 +2 5 6 +3 2 5 +4 5 8 +0 7 8 +4 4 7 +4 6 9 +2 8 9 + +vertices +10 + +patches + +knotvectors +2 +1 2 0 0 1 1 +1 2 0 0 1 1 + +dimension +2 + +controlpoints_cartesian +-1 -1 1 + 0 -1 1 +-1 0 1 + 0 0 1 + +knotvectors +2 +1 2 0 0 1 1 +1 2 0 0 1 1 + +dimension +2 + +controlpoints_cartesian +0 -1 1 +1 -1 1 +0 0 1 +1 0 1 + + +knotvectors +2 +1 2 0 0 1 1 +1 2 0 0 1 1 + +dimension +2 + +controlpoints_cartesian +-1 0 1 + 0 0 1 +-1 1 1 + 0 1 1 + +knotvectors +2 +1 2 0 0 1 1 +1 2 0 0 1 1 + +dimension +2 + +controlpoints_cartesian +0 0 1 +1 0 1 +0 1 1 +1 1 1 diff --git a/eltrans.cpp b/eltrans.cpp new file mode 100644 index 0000000..13b2bd5 --- /dev/null +++ b/eltrans.cpp @@ -0,0 +1,773 @@ +// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced +// at the Lawrence Livermore National Laboratory. All Rights reserved. See files +// LICENSE and NOTICE for details. LLNL-CODE-806117. +// +// This file is part of the MFEM library. For more information and source code +// availability visit https://mfem.org. +// +// MFEM is free software; you can redistribute it and/or modify it under the +// terms of the BSD-3 license. We welcome feedback and contributions, see file +// CONTRIBUTING.md for details. + +#include "../mesh/mesh_headers.hpp" +#include "fem.hpp" +#include "eltrans/eltrans_basis.hpp" + +#include + +namespace mfem +{ + +ElementTransformation::ElementTransformation() + : IntPoint(static_cast(NULL)), + EvalState(0), + geom(Geometry::INVALID), + Attribute(-1), + ElementNo(-1), + mesh(nullptr) +{ } + +real_t ElementTransformation::EvalWeight() +{ + MFEM_ASSERT((EvalState & WEIGHT_MASK) == 0, ""); + Jacobian(); + EvalState |= WEIGHT_MASK; + return (Wght = (dFdx.Width() == 0) ? 1.0 : dFdx.Weight()); +} + +const DenseMatrix &ElementTransformation::EvalAdjugateJ() +{ + MFEM_ASSERT((EvalState & ADJUGATE_MASK) == 0, ""); + Jacobian(); + adjJ.SetSize(dFdx.Width(), dFdx.Height()); + if (dFdx.Width() > 0) { CalcAdjugate(dFdx, adjJ); } + EvalState |= ADJUGATE_MASK; + return adjJ; +} + +const DenseMatrix &ElementTransformation::EvalTransAdjugateJ() +{ + MFEM_ASSERT((EvalState & TRANS_ADJUGATE_MASK) == 0, ""); + Jacobian(); + adjJT.SetSize(dFdx.Height(), dFdx.Width()); + if (dFdx.Width() == dFdx.Height()) { CalcAdjugateTranspose(dFdx, adjJT); } + else { AdjugateJacobian(); adjJT.Transpose(adjJ); } + EvalState |= TRANS_ADJUGATE_MASK; + return adjJT; +} + +const DenseMatrix &ElementTransformation::EvalInverseJ() +{ + // TODO: compute as invJ = / adjJ/Weight, if J is square, + // \ adjJ/Weight^2, otherwise. + MFEM_ASSERT((EvalState & INVERSE_MASK) == 0, ""); + Jacobian(); + invJ.SetSize(dFdx.Width(), dFdx.Height()); + if (dFdx.Width() > 0) { CalcInverse(dFdx, invJ); } + EvalState |= INVERSE_MASK; + return invJ; +} + +const DenseMatrix &ElementTransformation::EvalMetric() +{ + MFEM_ASSERT((EvalState & METRIC_MASK) == 0, ""); + + DenseMatrix JP(dFdx.Height()); + Geometries.JacToPerfJac(geom, Jacobian(), JP); + + Gij.SetSize(dFdx.Height(), dFdx.Height()); + MultAtB(JP, JP, Gij); + + EvalState |= METRIC_MASK; + return Gij; +} + +int InverseElementTransformation::FindClosestPhysPoint( + const Vector& pt, const IntegrationRule &ir) +{ + MFEM_VERIFY(T != NULL, "invalid ElementTransformation"); + MFEM_VERIFY(pt.Size() == T->GetSpaceDim(), "invalid point"); + + DenseMatrix physPts; + T->Transform(ir, physPts); + + // Initialize distance and index of closest point + int minIndex = -1; + real_t minDist = std::numeric_limits::max(); + + // Check all integration points in ir + const int npts = ir.GetNPoints(); + for (int i = 0; i < npts; ++i) + { + real_t dist = pt.DistanceTo(physPts.GetColumn(i)); + if (dist < minDist) + { + minDist = dist; + minIndex = i; + } + } + return minIndex; +} + +int InverseElementTransformation::FindClosestRefPoint( + const Vector& pt, const IntegrationRule &ir) +{ + MFEM_VERIFY(T != NULL, "invalid ElementTransformation"); + MFEM_VERIFY(pt.Size() == T->GetSpaceDim(), "invalid point"); + + // Initialize distance and index of closest point + int minIndex = -1; + real_t minDist = std::numeric_limits::max(); + + // Check all integration points in ir using the local metric at each point + // induced by the transformation. + Vector dp(T->GetSpaceDim()), dr(T->GetDimension()); + const int npts = ir.GetNPoints(); + for (int i = 0; i < npts; ++i) + { + const IntegrationPoint &ip = ir.IntPoint(i); + T->Transform(ip, dp); + dp -= pt; + T->SetIntPoint(&ip); + T->InverseJacobian().Mult(dp, dr); + real_t dist = dr.Norml2(); + // double dist = dr.Normlinf(); + if (dist < minDist) + { + minDist = dist; + minIndex = i; + } + } + return minIndex; +} + +void InverseElementTransformation::NewtonPrint(int mode, real_t val) +{ + std::ostream &os = mfem::out; + + // separator: + switch (mode%3) + { + case 0: os << ", "; break; + case 1: os << "Newton: "; break; + case 2: os << " "; break; + // "Newton: iter = xx, " + } + switch ((mode/3)%4) + { + case 0: os << "iter = " << std::setw(2) << int(val); break; + case 1: os << "delta_ref = " << std::setw(11) << val; break; + case 2: os << " err_phys = " << std::setw(11) << val; break; + case 3: break; + } + // ending: + switch ((mode/12)%4) + { + case 0: break; + case 1: os << '\n'; break; + case 2: os << " (converged)\n"; break; + case 3: os << " (actual)\n"; break; + } +} + +void InverseElementTransformation::NewtonPrintPoint(const char *prefix, + const Vector &pt, + const char *suffix) +{ + std::ostream &os = mfem::out; + + os << prefix << " = ("; + for (int j = 0; j < pt.Size(); j++) + { + os << (j > 0 ? ", " : "") << pt(j); + } + os << ')' << suffix; +} + +int InverseElementTransformation::NewtonSolve(const Vector &pt, + IntegrationPoint &ip) +{ + MFEM_ASSERT(pt.Size() == T->GetSpaceDim(), "invalid point"); + + const real_t phys_tol = phys_rtol*pt.Normlinf(); + + const int geom = T->GetGeometryType(); + const int dim = T->GetDimension(); + const int sdim = T->GetSpaceDim(); + IntegrationPoint xip, prev_xip; + real_t xd[3], yd[3], dxd[3], dxpd[3], dx_norm = -1.0, err_phys, + real_dx_norm = -1.0; + Vector x(xd, dim), y(yd, sdim), dx(dxd, dim), dx_prev(dxpd, dim); + bool hit_bdr = false, prev_hit_bdr = false; + + // Use ip0 as initial guess: + xip = ip0; + xip.Get(xd, dim); // xip -> x + if (print_level >= 3) + { + NewtonPrint(1, 0.); // iter 0 + NewtonPrintPoint(", ref_pt", x, "\n"); + } + + for (int it = 0; true; ) + { + // Remarks: + // If f(x) := 1/2 |pt-F(x)|^2, then grad(f)(x) = -J^t(x) [pt-F(x)]. + // Linearize F(y) at y=x: F(y) ~ L[x](y) := F(x) + J(x) [y-x]. + // Newton iteration for F(y)=b is given by L[x_old](x_new) = b, i.e. + // F(x_old) + J(x_old) [x_new-x_old] = b. + // + // To minimize: 1/2 |F(y)-b|^2, subject to: l(y) >= 0, we may consider the + // iteration: minimize: |L[x_old](x_new)-b|^2, subject to l(x_new) >= 0, + // i.e. minimize: |F(x_old) + J(x_old) [x_new-x_old] - b|^2. + + // This method uses: + // Newton iteration: x := x + J(x)^{-1} [pt-F(x)] + // or when dim != sdim: x := x + [J^t.J]^{-1}.J^t [pt-F(x)] + + // Compute the physical coordinates of the current point: + T->Transform(xip, y); + if (print_level >= 3) + { + NewtonPrint(11, 0.); // continuation line + NewtonPrintPoint("approx_pt", y, ", "); + NewtonPrintPoint("exact_pt", pt, "\n"); + } + subtract(pt, y, y); // y = pt-y + + // Check for convergence in physical coordinates: + err_phys = y.Normlinf(); + if (err_phys < phys_tol) + { + if (print_level >= 1) + { + NewtonPrint(1, (real_t)it); + NewtonPrint(3, dx_norm); + NewtonPrint(30, err_phys); + } + ip = xip; + if (solver_type != Newton) { return Inside; } + return Geometry::CheckPoint(geom, ip, ip_tol) ? Inside : Outside; + } + if (print_level >= 1) + { + if (it == 0 || print_level >= 2) + { + NewtonPrint(1, (real_t)it); + NewtonPrint(3, dx_norm); + NewtonPrint(18, err_phys); + } + } + + if (hit_bdr) + { + xip.Get(xd, dim); // xip -> x + if (prev_hit_bdr || it == max_iter || print_level >= 2) + { + prev_xip.Get(dxd, dim); // prev_xip -> dx + subtract(x, dx, dx); // dx = xip - prev_xip + real_dx_norm = dx.Normlinf(); + if (print_level >= 2) + { + NewtonPrint(41, real_dx_norm); + } + if (prev_hit_bdr && real_dx_norm < ref_tol) + { + if (print_level >= 0) + { + if (print_level <= 1) + { + NewtonPrint(1, (real_t)it); + NewtonPrint(3, dx_norm); + NewtonPrint(18, err_phys); + NewtonPrint(41, real_dx_norm); + } + mfem::out << "Newton: *** stuck on boundary!\n"; + } + return Outside; + } + } + } + + if (it == max_iter) { break; } + + // Perform a Newton step: + T->SetIntPoint(&xip); + T->InverseJacobian().Mult(y, dx); + x += dx; + it++; + if (solver_type != Newton) + { + prev_xip = xip; + prev_hit_bdr = hit_bdr; + } + xip.Set(xd, dim); // x -> xip + + // Perform projection based on solver_type: + switch (solver_type) + { + case Newton: break; + case NewtonSegmentProject: + hit_bdr = !Geometry::ProjectPoint(geom, prev_xip, xip); break; + case NewtonElementProject: + hit_bdr = !Geometry::ProjectPoint(geom, xip); break; + default: MFEM_ABORT("invalid solver type"); + } + if (print_level >= 3) + { + NewtonPrint(1, real_t(it)); + xip.Get(xd, dim); // xip -> x + NewtonPrintPoint(", ref_pt", x, "\n"); + } + + // Check for convergence in reference coordinates: + dx_norm = dx.Normlinf(); + if (dx_norm < ref_tol) + { + if (print_level >= 1) + { + NewtonPrint(1, (real_t)it); + NewtonPrint(27, dx_norm); + } + ip = xip; + if (solver_type != Newton) { return Inside; } + return Geometry::CheckPoint(geom, ip, ip_tol) ? Inside : Outside; + } + } + if (print_level >= 0) + { + if (print_level <= 1) + { + NewtonPrint(1, (real_t)max_iter); + NewtonPrint(3, dx_norm); + NewtonPrint(18, err_phys); + if (hit_bdr) { NewtonPrint(41, real_dx_norm); } + } + mfem::out << "Newton: *** iteration did not converge!\n"; + } + ip = xip; + return Unknown; +} + +int InverseElementTransformation::Transform(const Vector &pt, + IntegrationPoint &ip) +{ + MFEM_VERIFY(T != NULL, "invalid ElementTransformation"); + + // Select initial guess ... + switch (init_guess_type) + { + case Center: + ip0 = Geometries.GetCenter(T->GetGeometryType()); + break; + + case ClosestPhysNode: + case ClosestRefNode: + { + const int order = qpts_order >= 0 + ? qpts_order + : std::max(T->Order() + rel_qpts_order, 0); + if (order == 0) + { + ip0 = Geometries.GetCenter(T->GetGeometryType()); + } + else + { + RefinedGeometry &RefG = *refiner.Refine(T->GetGeometryType(), order); + int closest_idx = (init_guess_type == ClosestPhysNode) ? + FindClosestPhysPoint(pt, RefG.RefPts) : + FindClosestRefPoint(pt, RefG.RefPts); + ip0 = RefG.RefPts.IntPoint(closest_idx); + } + break; + } + case EdgeScan: + { + const int order = qpts_order >= 0 + ? qpts_order + : std::max(T->Order() + rel_qpts_order, 0); + if (order == 0) + { + ip0 = Geometries.GetCenter(T->GetGeometryType()); + } + else + { + auto &ir = *refiner.EdgeScan(T->GetGeometryType(), order + 1); + int res = Outside; + int npts = ir.GetNPoints(); + // will return Inside if any test point reports Inside, Outside if + // all points report Outside, else Unknown + for (int i = 0; i < npts; ++i) + { + ip0 = ir.IntPoint(i); + int tmp_res = NewtonSolve(pt, ip); + switch (tmp_res) + { + case Inside: + return Inside; + case Outside: + break; + case Unknown: + res = Unknown; + break; + } + } + return res; + } + break; + } + case GivenPoint: + break; + + default: + MFEM_ABORT("invalid initial guess type"); + } + + // Call the solver ... + return NewtonSolve(pt, ip); +} + + +void IsoparametricTransformation::SetIdentityTransformation( + Geometry::Type GeomType) +{ + switch (GeomType) + { + case Geometry::POINT : FElem = &PointFE; break; + case Geometry::SEGMENT : FElem = &SegmentFE; break; + case Geometry::TRIANGLE : FElem = &TriangleFE; break; + case Geometry::SQUARE : FElem = &QuadrilateralFE; break; + case Geometry::TETRAHEDRON : FElem = &TetrahedronFE; break; + case Geometry::CUBE : FElem = &HexahedronFE; break; + case Geometry::PRISM : FElem = &WedgeFE; break; + case Geometry::PYRAMID : FElem = &PyramidFE; break; + default: + MFEM_ABORT("unknown Geometry::Type!"); + } + int dim = FElem->GetDim(); + int dof = FElem->GetDof(); + const IntegrationRule &nodes = FElem->GetNodes(); + PointMat.SetSize(dim, dof); + for (int j = 0; j < dof; j++) + { + nodes.IntPoint(j).Get(&PointMat(0,j), dim); + } + geom = GeomType; +} + +const DenseMatrix &IsoparametricTransformation::EvalJacobian() +{ + MFEM_ASSERT((EvalState & JACOBIAN_MASK) == 0, ""); + + dshape.SetSize(FElem->GetDof(), FElem->GetDim()); + dFdx.SetSize(PointMat.Height(), dshape.Width()); + if (dshape.Width() > 0) + { + FElem->CalcDShape(*IntPoint, dshape); + Mult(PointMat, dshape, dFdx); + } + EvalState |= JACOBIAN_MASK; + + return dFdx; +} + +const DenseMatrix &IsoparametricTransformation::EvalHessian() +{ + MFEM_ASSERT((EvalState & HESSIAN_MASK) == 0, ""); + + int Dim = FElem->GetDim(); + d2shape.SetSize(FElem->GetDof(), (Dim*(Dim+1))/2); + d2Fdx2.SetSize(PointMat.Height(), d2shape.Width()); + if (d2shape.Width() > 0) + { + FElem->CalcHessian(*IntPoint, d2shape); + Mult(PointMat, d2shape, d2Fdx2); + } + EvalState |= HESSIAN_MASK; + + return d2Fdx2; +} + +int IsoparametricTransformation::OrderJ() const +{ + switch (FElem->Space()) + { + case FunctionSpace::Pk: + return (FElem->GetOrder()-1); + case FunctionSpace::Qk: + return (FElem->GetOrder()); + case FunctionSpace::Uk: + return (FElem->GetOrder()); + default: + MFEM_ABORT("unsupported finite element"); + } + return 0; +} + +int IsoparametricTransformation::OrderW() const +{ + switch (FElem->Space()) + { + case FunctionSpace::Pk: + return (FElem->GetOrder() - 1) * FElem->GetDim(); + case FunctionSpace::Qk: + return (FElem->GetOrder() * FElem->GetDim() - 1); + case FunctionSpace::Uk: + return (FElem->GetOrder() * FElem->GetDim() - 1); + default: + MFEM_ABORT("unsupported finite element"); + } + return 0; +} + +int IsoparametricTransformation::OrderGrad(const FiniteElement *fe) const +{ + if (FElem->Space() == fe->Space()) + { + int k = FElem->GetOrder(); + int d = FElem->GetDim(); + int l = fe->GetOrder(); + switch (fe->Space()) + { + case FunctionSpace::Pk: + return ((k-1)*(d-1)+(l-1)); + case FunctionSpace::Qk: + return (k*(d-1)+(l-1)); + case FunctionSpace::Uk: + return (k*(d-1)+(l-1)); + default: + MFEM_ABORT("unsupported finite element"); + } + } + MFEM_ABORT("incompatible finite elements"); + return 0; +} + +void IsoparametricTransformation::Transform (const IntegrationPoint &ip, + Vector &trans) +{ + MFEM_ASSERT(FElem != nullptr, "Must provide a valid FiniteElement object!"); + shape.SetSize(FElem->GetDof()); + trans.SetSize(PointMat.Height()); + + FElem -> CalcShape(ip, shape); + PointMat.Mult(shape, trans); +} + +void IsoparametricTransformation::Transform (const IntegrationRule &ir, + DenseMatrix &tr) +{ + int dof, n, dim, i, j, k; + + dim = PointMat.Height(); + dof = FElem->GetDof(); + n = ir.GetNPoints(); + + shape.SetSize(dof); + tr.SetSize(dim, n); + + for (j = 0; j < n; j++) + { + FElem -> CalcShape (ir.IntPoint(j), shape); + for (i = 0; i < dim; i++) + { + tr(i, j) = 0.0; + for (k = 0; k < dof; k++) + { + tr(i, j) += PointMat(i, k) * shape(k); + } + } + } +} + +void IsoparametricTransformation::Transform (const DenseMatrix &matrix, + DenseMatrix &result) +{ + MFEM_ASSERT(matrix.Height() == GetDimension(), "invalid input"); + result.SetSize(PointMat.Height(), matrix.Width()); + + IntegrationPoint ip; + Vector col; + + for (int j = 0; j < matrix.Width(); j++) + { + ip.Set(matrix.GetColumn(j), matrix.Height()); + + result.GetColumnReference(j, col); + Transform(ip, col); + } +} + +void IntegrationPointTransformation::Transform (const IntegrationPoint &ip1, + IntegrationPoint &ip2) +{ + real_t vec[3]; + Vector v (vec, Transf.GetPointMat().Height()); + + Transf.Transform (ip1, v); + ip2.Set(vec, v.Size()); +} + +void IntegrationPointTransformation::Transform (const IntegrationRule &ir1, + IntegrationRule &ir2) +{ + int i, n; + + n = ir1.GetNPoints(); + for (i = 0; i < n; i++) + { + Transform (ir1.IntPoint(i), ir2.IntPoint(i)); + } +} + +void FaceElementTransformations::SetIntPoint(const IntegrationPoint *face_ip) +{ + IsoparametricTransformation::SetIntPoint(face_ip); + + if (mask & 4) + { + Loc1.Transform(*face_ip, eip1); + if (Elem1) + { + Elem1->SetIntPoint(&eip1); + } + } + if (mask & 8) + { + Loc2.Transform(*face_ip, eip2); + if (Elem2) + { + Elem2->SetIntPoint(&eip2); + } + } +} + +ElementTransformation & +FaceElementTransformations::GetElement1Transformation() +{ + MFEM_VERIFY(mask & HAVE_ELEM1 && Elem1 != NULL, "The ElementTransformation " + "for the element has not been configured for side 1."); + return *Elem1; +} + +ElementTransformation & +FaceElementTransformations::GetElement2Transformation() +{ + MFEM_VERIFY(mask & HAVE_ELEM2 && Elem2 != NULL, "The ElementTransformation " + "for the element has not been configured for side 2."); + return *Elem2; +} + +IntegrationPointTransformation & +FaceElementTransformations::GetIntPoint1Transformation() +{ + MFEM_VERIFY(mask & HAVE_LOC1, "The IntegrationPointTransformation " + "for the element has not been configured for side 1."); + return Loc1; +} + +IntegrationPointTransformation & +FaceElementTransformations::GetIntPoint2Transformation() +{ + MFEM_VERIFY(mask & HAVE_LOC2, "The IntegrationPointTransformation " + "for the element has not been configured for side 2."); + return Loc2; +} + +void FaceElementTransformations::Transform(const IntegrationPoint &ip, + Vector &trans) +{ + MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation " + "for the face has not been configured."); + IsoparametricTransformation::Transform(ip, trans); +} + +void FaceElementTransformations::Transform(const IntegrationRule &ir, + DenseMatrix &tr) +{ + MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation " + "for the face has not been configured."); + IsoparametricTransformation::Transform(ir, tr); +} + +void FaceElementTransformations::Transform(const DenseMatrix &matrix, + DenseMatrix &result) +{ + MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation " + "for the face has not been configured."); + IsoparametricTransformation::Transform(matrix, result); +} + +real_t FaceElementTransformations::CheckConsistency(int print_level, + std::ostream &os) +{ + // Check that the face vertices are mapped to the same physical location + // when using the following three transformations: + // - the face transformation, *this + // - Loc1 + Elem1 + // - Loc2 + Elem2, if present. + + const bool have_face = (mask & 16); + const bool have_el1 = (mask & 1) && (mask & 4); + const bool have_el2 = (mask & 2) && (mask & 8) && (Elem2No >= 0); + if (int(have_face) + int(have_el1) + int(have_el2) < 2) + { + // need at least two different transformations to perform a check + return 0.0; + } + + const IntegrationRule &v_ir = *Geometries.GetVertices(GetGeometryType()); + + real_t max_dist = 0.0; + Vector dist(v_ir.GetNPoints()); + DenseMatrix coords_base, coords_el; + IntegrationRule v_eir(v_ir.GetNPoints()); + if (have_face) + { + Transform(v_ir, coords_base); + if (print_level > 0) + { + os << "\nface vertex coordinates (from face transform):\n" + << "----------------------------------------------\n"; + coords_base.PrintT(os, coords_base.Height()); + } + } + if (have_el1) + { + Loc1.Transform(v_ir, v_eir); + Elem1->Transform(v_eir, coords_el); + if (print_level > 0) + { + os << "\nface vertex coordinates (from element 1 transform):\n" + << "---------------------------------------------------\n"; + coords_el.PrintT(os, coords_el.Height()); + } + if (have_face) + { + coords_el -= coords_base; + coords_el.Norm2(dist); + max_dist = std::max(max_dist, dist.Normlinf()); + } + else + { + coords_base = coords_el; + } + } + if (have_el2) + { + Loc2.Transform(v_ir, v_eir); + Elem2->Transform(v_eir, coords_el); + if (print_level > 0) + { + os << "\nface vertex coordinates (from element 2 transform):\n" + << "---------------------------------------------------\n"; + coords_el.PrintT(os, coords_el.Height()); + } + coords_el -= coords_base; + coords_el.Norm2(dist); + max_dist = std::max(max_dist, dist.Normlinf()); + } + + return max_dist; +} +} diff --git a/eltrans.howto b/eltrans.howto new file mode 100644 index 0000000..4cc536f --- /dev/null +++ b/eltrans.howto @@ -0,0 +1,9 @@ +copy eltrans to mfem/fem +ans recompile + + +cp eltrans.* fem/fem + +cd build/mfem; make + +cd ../rbvms; make diff --git a/eltrans.hpp b/eltrans.hpp new file mode 100644 index 0000000..10c15e1 --- /dev/null +++ b/eltrans.hpp @@ -0,0 +1,929 @@ +// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced +// at the Lawrence Livermore National Laboratory. All Rights reserved. See files +// LICENSE and NOTICE for details. LLNL-CODE-806117. +// +// This file is part of the MFEM library. For more information and source code +// availability visit https://mfem.org. +// +// MFEM is free software; you can redistribute it and/or modify it under the +// terms of the BSD-3 license. We welcome feedback and contributions, see file +// CONTRIBUTING.md for details. + +#ifndef MFEM_ELEMENTTRANSFORM +#define MFEM_ELEMENTTRANSFORM + +#include "../config/config.hpp" +#include "../linalg/linalg.hpp" +#include "intrules.hpp" +#include "fe.hpp" + +#include "kernel_dispatch.hpp" + +namespace mfem +{ + +class GridFunction; + +class ElementTransformation +{ +protected: + const IntegrationPoint *IntPoint; + DenseMatrix dFdx, adjJ, invJ, Gij; + DenseMatrix d2Fdx2, adjJT; + real_t Wght; + int EvalState; + enum StateMasks + { + JACOBIAN_MASK = 1, + WEIGHT_MASK = 2, + ADJUGATE_MASK = 4, + INVERSE_MASK = 8, + HESSIAN_MASK = 16, + TRANS_ADJUGATE_MASK = 32, + METRIC_MASK = 64 + }; + Geometry::Type geom; + + /** @brief Evaluate the Jacobian of the transformation at the IntPoint and + store it in dFdx. */ + virtual const DenseMatrix &EvalJacobian() = 0; + + /** @brief Evaluate the Hessian of the transformation at the IntPoint and + store it in d2Fdx2. */ + virtual const DenseMatrix &EvalHessian() = 0; + + real_t EvalWeight(); + const DenseMatrix &EvalAdjugateJ(); + const DenseMatrix &EvalTransAdjugateJ(); + const DenseMatrix &EvalInverseJ(); + const DenseMatrix &EvalMetric(); + + /// @name Tolerance used for point comparisons + ///@{ +#ifdef MFEM_USE_DOUBLE + static constexpr real_t tol_0 = 1e-15; +#elif defined(MFEM_USE_SINGLE) + static constexpr real_t tol_0 = 1e-7; +#endif + ///@} + +public: + + /** This enumeration declares the values stored in + ElementTransformation::ElementType and indicates which group of objects + the index stored in ElementTransformation::ElementNo refers: + + | ElementType | Range of ElementNo + +-------------+------------------------- + | ELEMENT | [0, Mesh::GetNE() ) + | BDR_ELEMENT | [0, Mesh::GetNBE() ) + | EDGE | [0, Mesh::GetNEdges() ) + | FACE | [0, Mesh::GetNFaces() ) + | BDR_FACE | [0, Mesh::GetNBE() ) + */ + enum + { + ELEMENT = 1, + BDR_ELEMENT = 2, + EDGE = 3, + FACE = 4, + BDR_FACE = 5 + }; + + int Attribute, ElementNo, ElementType; + + /// The Mesh object containing the element. + /** If the element transformation belongs to a mesh, this will point to the + containing Mesh object. ElementNo will be the number of the element in + this Mesh. This will be NULL if the element does not belong to a mesh. */ + const Mesh *mesh; + + ElementTransformation(); + + /** @brief Force the reevaluation of the Jacobian in the next call. */ + void Reset() { EvalState = 0; } + + /** @brief Set the integration point @a ip that weights and Jacobians will + be evaluated at. */ + void SetIntPoint(const IntegrationPoint *ip) + { IntPoint = ip; EvalState = 0; } + + /** @brief Get a const reference to the currently set integration point. This + will return NULL if no integration point is set. */ + const IntegrationPoint &GetIntPoint() { return *IntPoint; } + + /** @brief Transform integration point from reference coordinates to + physical coordinates and store them in the vector. */ + virtual void Transform(const IntegrationPoint &, Vector &) = 0; + + /** @brief Transform all the integration points from the integration rule + from reference coordinates to physical + coordinates and store them as column vectors in the matrix. */ + virtual void Transform(const IntegrationRule &, DenseMatrix &) = 0; + + /** @brief Transform all the integration points from the column vectors + of @a matrix from reference coordinates to physical + coordinates and store them as column vectors in @a result. */ + virtual void Transform(const DenseMatrix &matrix, DenseMatrix &result) = 0; + + /** @brief Return the Jacobian matrix of the transformation at the currently + set IntegrationPoint, using the method SetIntPoint(). */ + /** The dimensions of the Jacobian matrix are physical-space-dim by + reference-space-dim. The first column contains the x derivatives of the + transformation, the second -- the y derivatives, etc. */ + const DenseMatrix &Jacobian() + { return (EvalState & JACOBIAN_MASK) ? dFdx : EvalJacobian(); } + + + /** @brief Return the Hessian matrix of the transformation at the currently + set IntegrationPoint, using the method SetIntPoint(). */ + const DenseMatrix &Hessian() + { return (EvalState & HESSIAN_MASK) ? d2Fdx2 : EvalHessian(); } + + /** @brief Return the weight of the Jacobian matrix of the transformation + at the currently set IntegrationPoint. + The Weight evaluates to $ \sqrt{\lvert J^T J \rvert} $. */ + real_t Weight() { return (EvalState & WEIGHT_MASK) ? Wght : EvalWeight(); } + + /** @brief Return the adjugate of the Jacobian matrix of the transformation + at the currently set IntegrationPoint. */ + const DenseMatrix &AdjugateJacobian() + { return (EvalState & ADJUGATE_MASK) ? adjJ : EvalAdjugateJ(); } + + /** @brief Return the transpose of the adjugate of the Jacobian matrix of + the transformation at the currently set IntegrationPoint. */ + const DenseMatrix &TransposeAdjugateJacobian() + { return (EvalState & TRANS_ADJUGATE_MASK) ? adjJT : EvalTransAdjugateJ(); } + + /** @brief Return the inverse of the Jacobian matrix of the transformation + at the currently set IntegrationPoint. */ + const DenseMatrix &InverseJacobian() + { return (EvalState & INVERSE_MASK) ? invJ : EvalInverseJ(); } + + /** @brief Return the metric tensor of the transformation + at the currently set IntegrationPoint. */ + const DenseMatrix &Metric() + { return (EvalState & METRIC_MASK) ? Gij : EvalMetric(); } + + /// Return the order of the current element we are using for the transformation. + virtual int Order() const = 0; + + /// Return the order of the elements of the Jacobian of the transformation. + virtual int OrderJ() const = 0; + + /** @brief Return the order of the determinant of the Jacobian (weight) + of the transformation. */ + virtual int OrderW() const = 0; + + /// Return the order of $ adj(J)^T \nabla fi $ + virtual int OrderGrad(const FiniteElement *fe) const = 0; + + /// Return the Geometry::Type of the reference element. + Geometry::Type GetGeometryType() const { return geom; } + + /// Return the topological dimension of the reference element. + int GetDimension() const { return Geometry::Dimension[geom]; } + + /// Get the dimension of the target (physical) space. + /** We support 2D meshes embedded in 3D; in this case the function will + return "3". */ + virtual int GetSpaceDim() const = 0; + + /** @brief Transform a point @a pt from physical space to a point @a ip in + reference space and optionally can set a solver tolerance using @a phys_tol. */ + /** Attempt to find the IntegrationPoint that is transformed into the given + point in physical space. If the inversion fails a non-zero value is + returned. This method is not 100 percent reliable for non-linear + transformations. */ + virtual int TransformBack(const Vector &pt, IntegrationPoint &ip, + const real_t phys_tol = tol_0) = 0; + + virtual ~ElementTransformation() { } +}; + + +/// The inverse transformation of a given ElementTransformation. +class InverseElementTransformation +{ +public: + /// Algorithms for selecting an initial guess. + enum InitGuessType + { + Center = 0, ///< Use the center of the reference element. + ClosestPhysNode = 1, /**< + Use the point returned by FindClosestPhysPoint() from a reference-space + grid of type and size controlled by SetInitGuessPointsType() and + SetInitGuessRelOrder(), respectively. */ + ClosestRefNode = 2, /**< + Use the point returned by FindClosestRefPoint() from a reference-space + grid of type and size controlled by SetInitGuessPointsType() and + SetInitGuessRelOrder(), respectively. */ + GivenPoint = 3, ///< Use a specific point, set with SetInitialGuess(). + EdgeScan = + 4, /**< Performs full solves on multiple points along the r/s/t=0 edges + of the element. It is recommended that SetInitGuessRelOrder() is + chosen such that max(trans_order+order,0)+1 <= 4 with + SetInitGuessPointsType() as Quadrature1D::ClosedUniform. @see + GeometryRefiner::EdgeScan */ + }; + + /// Solution strategy. + enum SolverType + { + Newton = 0, /**< + Use Newton's algorithm, without restricting the reference-space points + (iterates) to the reference element. */ + NewtonSegmentProject = 1, /**< + Use Newton's algorithm, restricting the reference-space points to the + reference element by scaling back the Newton increments, i.e. + projecting new iterates, x_new, lying outside the element, to the + intersection of the line segment [x_old, x_new] with the boundary. */ + NewtonElementProject = 2, /**< + Use Newton's algorithm, restricting the reference-space points to the + reference element by projecting new iterates, x_new, lying outside the + element, to the point on the boundary closest (in reference-space) to + x_new. */ + }; + + /// Values returned by Transform(). + enum TransformResult + { + Inside = 0, ///< The point is inside the element + Outside = 1, ///< The point is _probably_ outside the element + Unknown = 2 ///< The algorithm failed to determine where the point is + }; + +protected: + // Pointer to the forward transformation. Not owned. + ElementTransformation *T; + + // Parameters of the inversion algorithms: + IntegrationPoint ip0; + int init_guess_type; // algorithm to use + GeometryRefiner refiner; // geometry refiner for initial guess + int qpts_order; // num_1D_qpts = rel_qpts_order + 1, or < 0 to use + // rel_qpts_order. + int rel_qpts_order; // num_1D_qpts = max(trans_order+rel_qpts_order,0)+1 + int solver_type; // solution strategy to use + int max_iter; // max. number of Newton iterations + real_t ref_tol; // reference space tolerance + real_t phys_rtol; // physical space tolerance (relative) + real_t ip_tol; // tolerance for checking if a point is inside the ref. elem. + int print_level; + + void NewtonPrint(int mode, real_t val); + void NewtonPrintPoint(const char *prefix, const Vector &pt, + const char *suffix); + int NewtonSolve(const Vector &pt, IntegrationPoint &ip); + +public: + /// Construct the InverseElementTransformation with default parameters. + /** Some practical considerations regarding the choice of initial guess type + and solver type: + 1. The combination of #Center and #NewtonSegmentProject provides the + fastest way to estimate if a point lies inside an element, assuming + that most queried elements are not very deformed, e.g. if most + elements are convex. + 2. [Default] The combination of #Center and #NewtonElementProject + provides a somewhat slower alternative to 1 with the benefit of being + more reliable in the case when the query point is inside the element + but potentially slower in the case when the query point is outside the + element. + 3. The combination of #ClosestPhysNode and #NewtonElementProject is + slower than 1 and 2 but much more reliable, especially in the case of + highly distorted elements which do not have very high aspect ratios. + 4. The combination of #ClosestRefNode and #NewtonElementProject should + generally be the most reliable, coming at a bit higher computational + cost than 3 while performing better on distorted meshes with elements + having high aspect ratios. + @note None of these choices provide a guarantee that if a query point is + inside the element then it will be found. The only guarantee is that if + the Transform() method returns #Inside then the point lies inside the + element up to one of the specified physical- or reference-space + tolerances. */ + InverseElementTransformation(ElementTransformation *Trans = NULL) + : T(Trans), + init_guess_type(Center), + refiner(Quadrature1D::OpenHalfUniform), + qpts_order(-1), + rel_qpts_order(-1), + solver_type(NewtonElementProject), + max_iter(16), +#ifdef MFEM_USE_DOUBLE + ref_tol(1e-15), + phys_rtol(4e-15), + ip_tol(1e-8), +#elif defined(MFEM_USE_SINGLE) + ref_tol(4e-7), + phys_rtol(1e-6), + ip_tol(1e-4), +#endif + print_level(-1) + { } + + virtual ~InverseElementTransformation() { } + + /// Set a new forward ElementTransformation, @a Trans. + void SetTransformation(ElementTransformation &Trans) { T = &Trans; } + + /** @brief Choose how the initial guesses for subsequent calls to Transform() + will be selected. */ + void SetInitialGuessType(InitGuessType itype) { init_guess_type = itype; } + + /** @brief Set the initial guess for subsequent calls to Transform(), + switching to the #GivenPoint #InitGuessType at the same time. */ + void SetInitialGuess(const IntegrationPoint &init_ip) + { ip0 = init_ip; SetInitialGuessType(GivenPoint); } + + /// Set the Quadrature1D type used for the `Closest*` and `EdgeScan` initial + /// guess types. + void SetInitGuessPointsType(int q_type) { refiner.SetType(q_type); } + + /// Set the relative order used for the `Closest*` initial guess types. + /** The number of points in each spatial direction is given by the formula + max(trans_order+order,0)+1, where trans_order is the order of the current + ElementTransformation. */ + void SetInitGuessRelOrder(int order) + { + qpts_order = -1; + rel_qpts_order = order; + } + + /** The number of points in each spatial direction is given by the formula + order+1. */ + void SetInitGuessOrder(int order) + { + qpts_order = order; + } + + /** @brief Specify which algorithm to use for solving the transformation + equation, i.e. when calling the Transform() method. */ + void SetSolverType(SolverType stype) { solver_type = stype; } + + /// Set the maximum number of iterations when solving for a reference point. + void SetMaxIter(int max_it) { max_iter = max_it; } + + /// Set the reference-space convergence tolerance. + void SetReferenceTol(real_t ref_sp_tol) { ref_tol = ref_sp_tol; } + + /// Set the relative physical-space convergence tolerance. + void SetPhysicalRelTol(real_t phys_rel_tol) { phys_rtol = phys_rel_tol; } + + /** @brief Set the tolerance used to determine if a point lies inside or + outside of the reference element. */ + /** This tolerance is used only with the pure #Newton solver. */ + void SetElementTol(real_t el_tol) { ip_tol = el_tol; } + + /// Set the desired print level, useful for debugging. + /** The valid options are: -1 - never print (default); 0 - print only errors; + 1 - print the first and last iterations; 2 - print every iteration; + and 3 - print every iteration including point coordinates. */ + void SetPrintLevel(int pr_level) { print_level = pr_level; } + + /** @brief Find the IntegrationPoint mapped closest to @a pt. */ + /** This function uses the given IntegrationRule, @a ir, maps its points to + physical space and finds the one that is closest to the point @a pt. + + @param pt The query point. + @param ir The IntegrationRule, i.e. the set of reference points to map + to physical space and check. + @return The index of the IntegrationPoint in @a ir whose mapped point is + closest to @a pt. + @see FindClosestRefPoint(). */ + int FindClosestPhysPoint(const Vector& pt, const IntegrationRule &ir); + + /** @brief Find the IntegrationPoint mapped closest to @a pt, using a norm + that approximates the (unknown) distance in reference coordinates. */ + /** @see FindClosestPhysPoint(). */ + int FindClosestRefPoint(const Vector& pt, const IntegrationRule &ir); + + /** @brief Given a point, @a pt, in physical space, find its reference + coordinates, @a ip. + + @returns A value of type #TransformResult. */ + virtual int Transform(const Vector &pt, IntegrationPoint &ip); +}; + +/** + * @brief Performs batch inverse element transforms. Currently only supports + * non-mixed meshes with SEGMENT, SQUARE, or CUBE geometries. Mixed + * element order meshes are projected onto an equivalent uniform order mesh. + */ +class BatchInverseElementTransformation +{ + // nodes grid function, not owned + const GridFunction *gf_ = nullptr; + // initial guess algorithm to use + InverseElementTransformation::InitGuessType init_guess_type = + InverseElementTransformation::ClosestPhysNode; + int qpts_order = -1; // num_1D_qpts = rel_qpts_order + 1, or < 0 to use + // rel_qpts_order. + // num_1D_qpts = max(trans_order+rel_qpts_order,0)+1 + int rel_qpts_order = 0; + // solution strategy to use + InverseElementTransformation::SolverType solver_type = + InverseElementTransformation::NewtonElementProject; + // basis type stored in points1d + int basis_type = BasisType::Invalid; + // initial guess points type. Quadrature1D::Invalid is used for match + // basis_type. + int guess_points_type = Quadrature1D::Invalid; + // max. number of Newton iterations + int max_iter = 16; + // internal element node locations cache + Vector node_pos; +#ifdef MFEM_USE_DOUBLE + // reference space tolerance + real_t ref_tol = 1e-15; + // physical space tolerance (relative) + real_t phys_rtol = 4e-15; +#else + // reference space tolerance + real_t ref_tol = 4e-7; + // physical space tolerance (relative) + real_t phys_rtol = 1e-6; +#endif + // not owned, location of tensor product basis nodes in reference space + const Array *points1d = nullptr; + +public: + /// Uninitialized BatchInverseElementTransformation. Users must call + /// UpdateNodes before Transform. + BatchInverseElementTransformation(); + /// + /// Constructs a BatchInverseElementTransformation given @a nodes representing + /// the mesh nodes. + /// + BatchInverseElementTransformation(const GridFunction &nodes, + MemoryType d_mt = MemoryType::DEFAULT); + /// + /// Constructs a BatchInverseElementTransformation for a given @a mesh. + /// mesh.GetNodes() must not be null. + /// + BatchInverseElementTransformation(const Mesh &mesh, + MemoryType d_mt = MemoryType::DEFAULT); + + ~BatchInverseElementTransformation(); + + /** @brief Choose how the initial guesses for subsequent calls to Transform() + will be selected. ClosestRefNode is currently not supported. */ + void SetInitialGuessType(InverseElementTransformation::InitGuessType itype) + { + MFEM_ASSERT(itype != InverseElementTransformation::ClosestRefNode, + "ClosestRefNode is currently not supported"); + init_guess_type = itype; + } + + /// Set the Quadrature1D type used for the `Closest*` and `EdgeScan` initial + /// guess types. + void SetInitGuessPointsType(int q_type) { guess_points_type = q_type; } + + /// Set the relative order used for the `Closest*` initial guess types. + /** The number of points in each spatial direction is given by the formula + max(trans_order+order,0)+1, where trans_order is the order of the + current ElementTransformation. */ + void SetInitGuessRelOrder(int order) + { + qpts_order = -1; + rel_qpts_order = order; + } + + /** The number of points in each spatial direction is given by the formula + order+1. */ + void SetInitGuessOrder(int order) { qpts_order = order; } + + /// @b Gets the basis type nodes are projected onto, or BasisType::Invalid if + /// uninitialized. + int GetBasisType() const { return basis_type; } + + /** @brief Specify which algorithm to use for solving the transformation + equation, i.e. when calling the Transform() method. NewtonSegmentProject + is currently not supported. */ + void SetSolverType(InverseElementTransformation::SolverType stype) + { + MFEM_ASSERT(stype != InverseElementTransformation::NewtonSegmentProject, + "NewtonSegmentProject is currently not supported"); + solver_type = stype; + } + + /// Set the maximum number of iterations when solving for a reference point. + void SetMaxIter(int max_it) { max_iter = max_it; } + + /// Set the reference-space convergence tolerance. + void SetReferenceTol(real_t ref_sp_tol) { ref_tol = ref_sp_tol; } + + /// Set the relative physical-space convergence tolerance. + void SetPhysicalRelTol(real_t phys_rel_tol) { phys_rtol = phys_rel_tol; } + + /** + * @brief Updates internal datastructures if @a nodes change. Some version + * of UpdateNodes must be called at least once before calls to Transform if + * nodes have changed. + */ + void UpdateNodes(const GridFunction &nodes, + MemoryType d_mt = MemoryType::DEFAULT); + /** + * @brief Updates internal datastructures if @a mesh nodes change. Some version + * of UpdateNodes must be called at least once before calls to Transform if + * mesh nodes have changed. mesh.GetNodes() must not be null. + */ + void UpdateNodes(const Mesh &mesh, MemoryType d_mt = MemoryType::DEFAULT); + + /** @brief Performs a batch request of a set of points belonging to the given + elements. + @a pts list of physical point coordinates ordered by + Ordering::Type::byNODES. + @a elems which element index to search for each corresponding point in + @a pts + @a types output search classification (@see + InverseElementTransformation::TransformResult). + @a refs result reference point coordinates ordered by + Ordering::Type::byNODES. If using InitGuessType::GivenPoint, this should + contain the initial guess for each point. + @a use_device hint for if device acceleration should be used. + Device acceleration is currently only implemented for meshes containing + only a single tensor product basis element type. + @a iters optional array storing how many iterations was spent on each + tested point + */ + void Transform(const Vector &pts, const Array &elems, Array &types, + Vector &refs, bool use_device = true, + Array *iters = nullptr) const; + + using ClosestPhysPointKernelType = void (*)(int, int, int, int, + const real_t *, const real_t *, + const int *, const real_t *, + const real_t *, real_t *); + + // specialization params: Geom, SDim, use_device + MFEM_REGISTER_KERNELS(FindClosestPhysPoint, ClosestPhysPointKernelType, + (int, int, bool)); + + using ClosestPhysDofKernelType = void (*)(int, int, int, + const real_t *, const real_t *, + const int *, const real_t *, + real_t *); + + // specialization params: Geom, SDim, use_device + MFEM_REGISTER_KERNELS(FindClosestPhysDof, ClosestPhysDofKernelType, + (int, int, bool)); + + using ClosestRefDofKernelType = void (*)(int, int, int, const real_t *, + const real_t *, const int *, + const real_t *, real_t *); + + // specialization params: Geom, SDim, use_device + MFEM_REGISTER_KERNELS(FindClosestRefDof, ClosestRefDofKernelType, + (int, int, bool)); + + using ClosestRefPointKernelType = void (*)(int, int, int, int, const real_t *, + const real_t *, const int *, + const real_t *, const real_t *, + real_t *); + + // specialization params: Geom, SDim, use_device + MFEM_REGISTER_KERNELS(FindClosestRefPoint, ClosestRefPointKernelType, + (int, int, bool)); + + using NewtonKernelType = void (*)(real_t, real_t, int, int, int, int, + const real_t *, const real_t *, + const int *, const real_t *, int *, int*, + real_t *); + + // specialization params: Geom, SDim, SolverType, use_device + MFEM_REGISTER_KERNELS(NewtonSolve, NewtonKernelType, + (int, int, InverseElementTransformation::SolverType, + bool)); + + using NewtonEdgeScanKernelType = void (*)(real_t, real_t, int, int, int, int, + const real_t *, const real_t *, + const int *, const real_t *, + const real_t *, int, int *, int *, + real_t *); + + // specialization params: Geom, SDim, SolverType, use_device + MFEM_REGISTER_KERNELS(NewtonEdgeScan, NewtonEdgeScanKernelType, + (int, int, InverseElementTransformation::SolverType, + bool)); + + struct Kernels { Kernels(); }; + + template + static void AddFindClosestSpecialization() + { + FindClosestPhysPoint::Specialization::Add(); + FindClosestRefPoint::Specialization::Add(); + FindClosestPhysPoint::Specialization::Add(); + FindClosestRefPoint::Specialization::Add(); + FindClosestPhysDof::Specialization::Add(); + FindClosestRefDof::Specialization::Add(); + FindClosestPhysDof::Specialization::Add(); + FindClosestRefDof::Specialization::Add(); + } + + template + static void AddNewtonSolveSpecialization() + { + NewtonSolve::Specialization::Add(); + NewtonEdgeScan::Specialization::Add(); + NewtonSolve::Specialization::Add(); + NewtonEdgeScan::Specialization::Add(); + } +}; + +/// A standard isoparametric element transformation +class IsoparametricTransformation : public ElementTransformation +{ +private: + DenseMatrix dshape, d2shape; + Vector shape; + + const FiniteElement *FElem; + DenseMatrix PointMat; // dim x dof + + /** @brief Evaluate the Jacobian of the transformation at the IntPoint and + store it in dFdx. */ + const DenseMatrix &EvalJacobian() override; + // Evaluate the Hessian of the transformation at the IntPoint and store it + // in d2Fdx2. + const DenseMatrix &EvalHessian() override; + +public: + IsoparametricTransformation() : FElem(NULL) {} + + /// Set the element that will be used to compute the transformations + void SetFE(const FiniteElement *FE) + { + MFEM_ASSERT(FE != NULL, "Must provide a valid FiniteElement object!"); + EvalState = (FE != FElem) ? 0 : EvalState; + FElem = FE; geom = FE->GetGeomType(); + } + + /// Get the current element used to compute the transformations + const FiniteElement* GetFE() const { return FElem; } + + /// @brief Set the underlying point matrix describing the transformation. + /** The dimensions of the matrix are space-dim x dof. The transformation is + defined as + $ x = F( \hat x ) = P \phi( \hat x ) $ + + where $ \hat x $ is the reference point, @a x is the corresponding + physical point, @a P is the point matrix, and $ \phi( \hat x ) $ is + the column-vector of all basis functions evaluated at $ \hat x $ . + The columns of @a P represent the control points in physical space + defining the transformation. */ + void SetPointMat(const DenseMatrix &pm) { PointMat = pm; EvalState = 0; } + + /// Return the stored point matrix. + const DenseMatrix &GetPointMat() const { return PointMat; } + + /// @brief Write access to the stored point matrix. Use with caution. + /** If the point matrix is altered using this member function the Reset + function should also be called to force the reevaluation of the + Jacobian, etc.. */ + DenseMatrix &GetPointMat() { return PointMat; } + + /// Set the FiniteElement Geometry for the reference elements being used. + void SetIdentityTransformation(Geometry::Type GeomType); + + /** @brief Transform integration point from reference coordinates to + physical coordinates and store them in the vector. */ + void Transform(const IntegrationPoint &, Vector &) override; + + /** @brief Transform all the integration points from the integration rule + from reference coordinates to physical + coordinates and store them as column vectors in the matrix. */ + void Transform(const IntegrationRule &, DenseMatrix &) override; + + /** @brief Transform all the integration points from the column vectors + of @a matrix from reference coordinates to physical + coordinates and store them as column vectors in @a result. */ + void Transform(const DenseMatrix &matrix, DenseMatrix &result) override; + + /// Return the order of the current element we are using for the transformation. + int Order() const override { return FElem->GetOrder(); } + + /// Return the order of the elements of the Jacobian of the transformation. + int OrderJ() const override; + + /** @brief Return the order of the determinant of the Jacobian (weight) + of the transformation. */ + int OrderW() const override; + + /// Return the order of $ adj(J)^T \nabla fi $ + int OrderGrad(const FiniteElement *fe) const override; + + int GetSpaceDim() const override { return PointMat.Height(); } + + /** @brief Transform a point @a pt from physical space to a point @a ip in + reference space and optionally can set a solver tolerance using @a phys_tol. */ + /** Attempt to find the IntegrationPoint that is transformed into the given + point in physical space. If the inversion fails a non-zero value is + returned. This method is not 100 percent reliable for non-linear + transformations. */ + int TransformBack (const Vector & v, IntegrationPoint & ip, + const real_t phys_rel_tol = tol_0) override + { + InverseElementTransformation inv_tr(this); + inv_tr.SetPhysicalRelTol(phys_rel_tol); + return inv_tr.Transform(v, ip); + } + + virtual ~IsoparametricTransformation() { } + + MFEM_DEPRECATED void FinalizeTransformation() {} +}; + +class IntegrationPointTransformation +{ +public: + IsoparametricTransformation Transf; + void Transform (const IntegrationPoint &, IntegrationPoint &); + void Transform (const IntegrationRule &, IntegrationRule &); +}; + +/** @brief A specialized ElementTransformation class representing a face and + its two neighboring elements. + + This class can be used as a container for the element transformation data + needed for integrating discontinuous fields on element interfaces in a + Discontinuous Galerkin (DG) context. + + The secondary purpose of this class is to enable the + GridFunction::GetValue function, and various related functions, to properly + evaluate fields with limited continuity on boundary elements. +*/ +class FaceElementTransformations : public IsoparametricTransformation +{ +private: + + // Bitwise OR of ConfigMasks + int mask; + + IntegrationPoint eip1, eip2; + +protected: // interface for Mesh to be able to configure this object. + + friend class Mesh; +#ifdef MFEM_USE_MPI + friend class ParMesh; +#endif + + /// Set the mask indicating which portions of the object have been setup + /** The argument @a m is a bitmask used in + Mesh::GetFaceElementTransformations to indicate which portions of the + FaceElementTransformations object have been configured. + + mask & 1: Elem1 is configured + mask & 2: Elem2 is configured + mask & 4: Loc1 is configured + mask & 8: Loc2 is configured + mask & 16: The Face transformation itself is configured + */ + void SetConfigurationMask(int m) { mask = m; } + +public: + + enum ConfigMasks + { + HAVE_ELEM1 = 1, ///< Element on side 1 is configured + HAVE_ELEM2 = 2, ///< Element on side 2 is configured + HAVE_LOC1 = 4, ///< Point transformation for side 1 is configured + HAVE_LOC2 = 8, ///< Point transformation for side 2 is configured + HAVE_FACE = 16 ///< Face transformation is configured + }; + + int Elem1No, Elem2No; + Geometry::Type &FaceGeom; ///< @deprecated Use GetGeometryType instead + ElementTransformation *Elem1, *Elem2; + ElementTransformation *Face; ///< @deprecated No longer necessary + IntegrationPointTransformation Loc1, Loc2; + + FaceElementTransformations() : FaceGeom(geom), Face(this) {} + + /** @brief Method to set the geometry type of the face. + + @note This method is designed to be used when + [Par]Mesh::GetFaceTransformation will not be called i.e. when the face + transformation will not be needed but the neighboring element + transformations will be. Using this method to override the GeometryType + should only be done with great care. + */ + void SetGeometryType(Geometry::Type g) { geom = g; } + + /** @brief Return the mask defining the configuration state. + + The mask value indicates which portions of FaceElementTransformations + object have been configured. + + mask & 1: Elem1 is configured + mask & 2: Elem2 is configured + mask & 4: Loc1 is configured + mask & 8: Loc2 is configured + mask & 16: The Face transformation itself is configured + */ + int GetConfigurationMask() const { return mask; } + + /** @brief Set the integration point in the Face and the two neighboring + elements, if present. + + The point @a face_ip must be in the reference coordinate system of the + face. + */ + void SetIntPoint(const IntegrationPoint *face_ip); + + /** @brief Set the integration point in the Face and the two neighboring + elements, if present. + + This is a more expressive member function name than SetIntPoint, which + in this special case, does the same thing. This function can be used for + greater code clarity. + */ + inline void SetAllIntPoints(const IntegrationPoint *face_ip) + { FaceElementTransformations::SetIntPoint(face_ip); } + + /** @brief Get a const reference to the integration point in neighboring + element 1 corresponding to the currently set integration point on the + face. + + This IntegrationPoint object will only contain up-to-date data if + SetIntPoint or SetAllIntPoints has been called with the latest + integration point for the face and the appropriate point transformation + has been configured. */ + const IntegrationPoint &GetElement1IntPoint() { return eip1; } + + /** @brief Get a const reference to the integration point in neighboring + element 2 corresponding to the currently set integration point on the + face. + + This IntegrationPoint object will only contain up-to-date data if + SetIntPoint or SetAllIntPoints has been called with the latest + integration point for the face and the appropriate point transformation + has been configured. */ + const IntegrationPoint &GetElement2IntPoint() { return eip2; } + + void Transform(const IntegrationPoint &, Vector &) override; + void Transform(const IntegrationRule &, DenseMatrix &) override; + void Transform(const DenseMatrix &matrix, DenseMatrix &result) override; + + ElementTransformation & GetElement1Transformation(); + ElementTransformation & GetElement2Transformation(); + IntegrationPointTransformation & GetIntPoint1Transformation(); + IntegrationPointTransformation & GetIntPoint2Transformation(); + + /** @brief Check for self-consistency: compares the result of mapping the + reference face vertices to physical coordinates using the three + transformations: face, element 1, and element 2. + + @param[in] print_level If set to a positive number, print the physical + coordinates of the face vertices computed through + all available transformations: face, element 1, + and/or element 2. + @param[in,out] out The output stream to use for printing. + + @returns A maximal distance between physical coordinates of face vertices + that should coincide. A successful check should return a small + number relative to the mesh extents. If less than 2 of the three + transformations are set, returns 0. + + @warning This check will generally fail on periodic boundary faces. + */ + real_t CheckConsistency(int print_level = 0, + std::ostream &out = mfem::out); +}; + +/** Elem1(Loc1(x)) = Face(x) = Elem2(Loc2(x)) + + + Physical Space + + *--------* ^ *--------* + Elem1No / / \ / \ / \ \ Elem2No + / / \ / \ / \ \ + / / n \ / \ / \ \ + *--------* ==> * ( ) * *--------* + \ \ / \ / \ / / + \ \ / \ / \ / / + \ \ / \ / \ / / + *--------* v *--------* + + ^ ^ + | ^ | + Elem1 | | | Elem2 + | | Face | + | + *--------* *--------* + / /| / /| + 1 *--------* | 1 *--------* 1 *--------* | + | | | Loc1 | | Loc2 | | | + | | * <----- | x | -----> | | * + | |/ | | | |/ + *--------* *--------* *--------* + 0 1 0 1 0 1 + + Reference Space +*/ + +} + +#endif diff --git a/fe_base.patch b/fe_base.patch new file mode 100644 index 0000000..81ef841 --- /dev/null +++ b/fe_base.patch @@ -0,0 +1,14 @@ +diff --git a/fem/fe/fe_base.cpp b/fem/fe/fe_base.cpp +index ab1977ddbe..c350cae54c 100644 +--- a/fem/fe/fe_base.cpp ++++ b/fem/fe/fe_base.cpp +@@ -111,7 +111,8 @@ void FiniteElement::GetFaceDofs(int face, int **dofs, int *ndofs) const + void FiniteElement::CalcHessian(const IntegrationPoint &ip, + DenseMatrix &h) const + { +- MFEM_ABORT("method is not overloaded"); ++ h = 0.0; ++ // MFEM_ABORT("method is not overloaded"); + } + + void FiniteElement::GetLocalInterpolation(ElementTransformation &Trans, diff --git a/supg/CMakeLists.txt b/supg/CMakeLists.txt index b99d705..b6c1af9 100644 --- a/supg/CMakeLists.txt +++ b/supg/CMakeLists.txt @@ -27,3 +27,13 @@ endif() add_executable(supg ${SOURCES} ${HEADERS}) target_link_libraries(supg PRIVATE ${LIBRARIES} util ${CMAKE_DL_LIBS}) install(TARGETS supg DESTINATION bin) + +# Small test utility to check tau variance under vertex-order permutations +add_executable(tau_variance tau_variance.cpp) +target_link_libraries(tau_variance PRIVATE ${LIBRARIES} util ${CMAKE_DL_LIBS}) +install(TARGETS tau_variance DESTINATION bin) + +# Alternate variant with direction sweeps and separate convective/diffusive checks +add_executable(tau_variance_alt tau_variance_alt.cpp) +target_link_libraries(tau_variance_alt PRIVATE ${LIBRARIES} util ${CMAKE_DL_LIBS}) +install(TARGETS tau_variance_alt DESTINATION bin) diff --git a/supg/integrator.cpp b/supg/integrator.cpp index 0692685..b07e734 100644 --- a/supg/integrator.cpp +++ b/supg/integrator.cpp @@ -9,6 +9,149 @@ using namespace mfem; +// This file is part of the RBVMS application. For more information and source +// code availability visit https://idoakkerman.github.io/ +// +// RBVMS is free software; you can redistribute it and/or modify it under the +// terms of the BSD-3 license. +//------------------------------------------------------------------------------ + +#include "integrator.hpp" + +using namespace mfem; + +// Define the integration rule based on FE order +const IntegrationRule &StabTauIntegrator::GetRule( + const FiniteElement &trial_fe, + const FiniteElement &test_fe, + ElementTransformation &Trans) +{ + int order = trial_fe.GetOrder() + test_fe.GetOrder(); + return IntRules.Get(trial_fe.GetGeomType(), order); +} + +// Define convective tau +real_t StabTauIntegrator::GetTau(real_t &k, Vector &a, DenseMatrix &Gij, + real_t Ch) +{ + // printf("C*h = %f\n", Ch); + Ch = 1/Ch; + double tau = 1e-10; + int dim = Gij.Width(); + for (int j = 0; j < dim; j++) + { + for (int i = 0; i < dim; i++) + { + tau += Gij(i,j)*a[i]*a[j] + Ch*Ch*k*k; + } + } + return 1.0/sqrt(tau); +} + + + +// Set the dimension of temporary variables +void StabTauIntegrator::SetDim(int d) +{ + if (dim != d ) + { + dim = d; + Gij.SetSize(dim); + a.SetSize(dim); + dphidx.SetSize(dim); + } +} + +// Compute the element nonlinear residual +void StabTauIntegrator::AssembleElementVector(const FiniteElement &el, + ElementTransformation &Trans, + const Vector &elfun, + Vector &elvect) +{ + real_t w,mu,f,phi,dphidx2,tau; + + SetDim(el.GetDim()); + int nd = el.GetDof(); + + elvect.SetSize(nd); + shape.SetSize(nd); + dshape.SetSize(nd,dim); + lshape.SetSize(nd); + const IntegrationRule *ir = NonlinearFormIntegrator::IntRule ? + NonlinearFormIntegrator::IntRule : &GetRule(el, el, Trans); + + elvect = 0.0; + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + Trans.SetIntPoint (&ip); + w = Trans.Weight() * ip.weight; + + //adv_cf->Eval(a, Trans, ip); + mu = mu_cf->Eval(Trans, ip); + + // Calculate shapes + el.CalcPhysShape(Trans, shape); + tau = shape*elfun; + el.CalcPhysDShape(Trans, dshape); + el.CalcPhysLaplacian(Trans, lshape); + for (int j = 0; j < nd; j++) + { + double val = mu * mu* lshape(j)*lshape(j)*tau; + for (int d = 0; d < dim; d++) + { + val -= mu*dshape(j,d) * dshape(j,d); + } + elvect[j] += val*w; + } + } + //elvect.Print(); +} + +// Compute the element jacobian +void StabTauIntegrator::AssembleElementGrad(const FiniteElement &el, + ElementTransformation &Trans, + const Vector &elfun, + DenseMatrix &elmat) +{ + + real_t w,mu,f,phi,dphidx2,tau; + + SetDim(el.GetDim()); + int nd = el.GetDof(); + double num = 0.0, den = 0.0; + + elmat.SetSize(nd); + shape.SetSize(nd); + lshape.SetSize(nd); + const IntegrationRule *ir = NonlinearFormIntegrator::IntRule ? + NonlinearFormIntegrator::IntRule : &GetRule(el, el, Trans); + + elmat = 0.0; + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + Trans.SetIntPoint (&ip); + w = Trans.Weight() * ip.weight; + + //adv_cf->Eval(a, Trans, ip); + mu = mu_cf->Eval(Trans, ip); + + // Calculate shapes + el.CalcPhysShape(Trans, shape); + el.CalcPhysLaplacian(Trans, lshape); + + // Laplacian + for (int j = 0; j < nd; j++) + { + for (int k = 0; k < nd; k++) + { + elmat(j,k) += w * mu * mu * lshape(j)*lshape(j)*shape(k); + } + } + } +} + // Define the integration rule based on FE order const IntegrationRule &StabConvDifIntegrator::GetRule( const FiniteElement &trial_fe, @@ -20,16 +163,18 @@ const IntegrationRule &StabConvDifIntegrator::GetRule( } // Define convective tau -real_t StabConvDifIntegrator::GetTau(real_t &k, Vector &a, DenseMatrix &Gij) +real_t StabConvDifIntegrator::GetTau(real_t &k, Vector &a, DenseMatrix &Gij, + real_t Ch) { - double CI = 1./12.0; + // printf("C*h = %f\n", Ch); + Ch = 1/Ch; double tau = 1e-10; int dim = Gij.Width(); for (int j = 0; j < dim; j++) { for (int i = 0; i < dim; i++) { - tau += Gij(i,j)*a[i]*a[j] + CI*CI*k*k*Gij(i,j)*Gij(i,j); + tau += Gij(i,j)*a[i]*a[j] + Ch*Ch*k*k; } } return 1.0/sqrt(tau); @@ -37,9 +182,9 @@ real_t StabConvDifIntegrator::GetTau(real_t &k, Vector &a, DenseMatrix &Gij) // Define convective kdc real_t StabConvDifIntegrator::GetKdc(Vector &a, - real_t &res, - Vector &dphidx, - DenseMatrix &Gij) + real_t &res, + Vector &dphidx, + DenseMatrix &Gij) { real_t h = 1.0/sqrt(Gij.Trace()/Gij.Width()); return kdc0*h*a.Norml2() + kdc1*h*fabs(res)/(dphidx.Norml2() + 1e-10); @@ -59,9 +204,9 @@ void StabConvDifIntegrator::SetDim(int d) // Compute the element nonlinear residual void StabConvDifIntegrator::AssembleElementVector(const FiniteElement &el, - ElementTransformation &Trans, - const Vector &elfun, - Vector &elvect) + ElementTransformation &Trans, + const Vector &elfun, + Vector &elvect) { real_t w,mu,f,phi,dphidx2,res,res_red; @@ -83,7 +228,16 @@ void StabConvDifIntegrator::AssembleElementVector(const FiniteElement &el, const IntegrationPoint &ip = ir->IntPoint(i); Trans.SetIntPoint (&ip); w = Trans.Weight() * ip.weight; - MultAtB(Trans.InverseJacobian(),Trans.InverseJacobian(),Gij); + + // Metric tensor (with perfect geometry) + const mfem::DenseMatrix &A = Geometries.GetGeomToPerfGeomJac(Trans.GetGeometryType()); + const mfem::DenseMatrix &invJ = Trans.InverseJacobian(); + mfem::DenseMatrix invJ_perf(invJ.Height(), invJ.Width()); + Mult(A, invJ, invJ_perf); + MultAtB(invJ_perf, invJ_perf, Gij); + + // Old metric tensor derivation: + // MultAtB(Trans.InverseJacobian(), Trans.InverseJacobian(), Gij); // Calculate shapes el.CalcPhysShape(Trans, shape); @@ -114,15 +268,16 @@ void StabConvDifIntegrator::AssembleElementVector(const FiniteElement &el, test.Add((int) type*mu, lshape); // Add Reaction stabilization term // Stabilized terms - elvect.Add(w*GetTau(mu, a, Gij)*res, test); + real_t Ch = inv_cf->Eval(Trans, ip); + elvect.Add(w*GetTau(mu, a, Gij, Ch)*res, test); } } // Compute the element jacobian void StabConvDifIntegrator::AssembleElementGrad(const FiniteElement &el, - ElementTransformation &Trans, - const Vector &elfun, - DenseMatrix &elmat) + ElementTransformation &Trans, + const Vector &elfun, + DenseMatrix &elmat) { real_t w,mu,f,phi,dphidx2,res; @@ -145,7 +300,16 @@ void StabConvDifIntegrator::AssembleElementGrad(const FiniteElement &el, const IntegrationPoint &ip = ir->IntPoint(i); Trans.SetIntPoint (&ip); w = Trans.Weight() * ip.weight; - MultAtB(Trans.InverseJacobian(),Trans.InverseJacobian(),Gij); + + // Metric tensor (with perfect geometry) + const mfem::DenseMatrix &A = Geometries.GetGeomToPerfGeomJac(Trans.GetGeometryType()); + const mfem::DenseMatrix &invJ = Trans.InverseJacobian(); + mfem::DenseMatrix invJ_perf(invJ.Height(), invJ.Width()); + Mult(A, invJ, invJ_perf); + MultAtB(invJ_perf, invJ_perf, Gij); + + // Old metric tensor derivation + // MultAtB(Trans.InverseJacobian(), Trans.InverseJacobian(), Gij); // Calculate shapes el.CalcPhysShape(Trans, shape); @@ -178,6 +342,7 @@ void StabConvDifIntegrator::AssembleElementGrad(const FiniteElement &el, dshape.Mult(a, trail); // Add Convection term trail.Add(-mu, lshape); // Add Diffusion term - AddMult_a_VWt(w*GetTau(mu, a, Gij), test, trail, elmat); + real_t Ch = inv_cf->Eval(Trans, ip); + AddMult_a_VWt(w*GetTau(mu, a, Gij, Ch), test, trail, elmat); } } diff --git a/supg/integrator.hpp b/supg/integrator.hpp index 9a3ac61..b7be2a2 100644 --- a/supg/integrator.hpp +++ b/supg/integrator.hpp @@ -12,12 +12,66 @@ #ifndef SUPG_INTEGRATOR_HPP #define SUPG_INTEGRATOR_HPP +#include "coefficients.hpp" #include "mfem.hpp" using namespace mfem; +/** This Class defines an integrator for +*/ +class StabTauIntegrator : public NonlinearFormIntegrator +{ +private: + // Physical parameters + VectorCoefficient *adv_cf; + Coefficient *mu_cf; + + InverseEstimateCoefficient *inv_cf; + + /// The stabilization parameter + int dim; + void SetDim(int dim); + + /// The stabilization parameter + real_t GetTau(real_t &k, Vector &a, DenseMatrix &Gij, real_t CI); + + /// Temporary variables + Vector a, dphidx, shape, lshape, trail, test; + DenseMatrix dshape, Gij; + +public: + /// Constructor + StabTauIntegrator(VectorCoefficient &a, + Coefficient &m, + InverseEstimateCoefficient &c) + : adv_cf(&a), mu_cf(&m), inv_cf(&c) + { + dim = -1; + }; -enum StabilizeType{ + /// Destructor + ~StabTauIntegrator() {}; + + /// Compute the element nonlinear residual + virtual void AssembleElementVector(const FiniteElement &el, + ElementTransformation &Tr, + const Vector &elfun, + Vector &elvect) override; + + /// Compute the element jacobian + virtual void AssembleElementGrad(const FiniteElement &el, + ElementTransformation &Tr, + const Vector &elfun, + DenseMatrix &elmat) override; + + /// Give element integration rule + static const IntegrationRule &GetRule(const FiniteElement &trial_fe, + const FiniteElement &test_fe, + ElementTransformation &Trans); +}; + +enum StabilizeType +{ GLS = -1, SUPG = 0, VMS = 1 @@ -37,13 +91,15 @@ class StabConvDifIntegrator : public NonlinearFormIntegrator Coefficient *mu_cf; Coefficient *force_cf; + InverseEstimateCoefficient *inv_cf; + /// The stabilization parameter int dim; void SetDim(int dim); /// The stabilization parameter StabilizeType type; - real_t GetTau(real_t &k, Vector &a, DenseMatrix &Gij); + real_t GetTau(real_t &k, Vector &a, DenseMatrix &Gij, real_t CI); /// The discontinuity capturing parameter real_t kdc0; // inconsistent part @@ -59,8 +115,9 @@ class StabConvDifIntegrator : public NonlinearFormIntegrator StabConvDifIntegrator(VectorCoefficient &a, Coefficient &m, Coefficient &f, + InverseEstimateCoefficient &c, real_t k0 = 0.0, - real_t k1 = 0.0) : adv_cf(&a), mu_cf(&m), force_cf(&f) + real_t k1 = 0.0) : adv_cf(&a), mu_cf(&m), force_cf(&f), inv_cf(&c) { type = StabilizeType::SUPG; kdc0 = k0; @@ -68,10 +125,10 @@ class StabConvDifIntegrator : public NonlinearFormIntegrator dim = -1; }; - void SetStabilization(StabilizeType t){ type = t; }; - void SetGLS(){ type = StabilizeType::GLS; }; - void SetSUPG(){ type = StabilizeType::SUPG; }; - void SetVMS(){ type = StabilizeType::VMS; }; + void SetStabilization(StabilizeType t) { type = t; }; + void SetGLS() { type = StabilizeType::GLS; }; + void SetSUPG() { type = StabilizeType::SUPG; }; + void SetVMS() { type = StabilizeType::VMS; }; /// Destructor ~StabConvDifIntegrator() {}; diff --git a/supg/supg.cpp b/supg/supg.cpp index 963d28f..ec93a0e 100644 --- a/supg/supg.cpp +++ b/supg/supg.cpp @@ -143,6 +143,7 @@ int main(int argc, char *argv[]) // Refine mesh { + // mesh.DegreeElevate(1); if (mesh.NURBSext && (strlen(ref_file) != 0)) { mesh.RefineNURBSFromFile(ref_file); @@ -223,103 +224,196 @@ int main(int argc, char *argv[]) } } } + /* + + + TAU PART + + */ - // Define the gridfunction and solution vector - ParGridFunction phi_gf(space); - LibCoefficient sol_phi(lib_file, "sol_phi"); - phi_gf.ProjectCoefficient(sol_phi); - Vector xp; - phi_gf.GetTrueDofs(xp); - - // Define the visualisation output - VisItDataCollection vdc("step", &pmesh); - vdc.SetPrefixPath(vis_dir); - vdc.RegisterField("phi", &phi_gf); - vdc.SetCycle(0); - vdc.Save(); - - // Define the physical parameters - LibVectorCoefficient adv(dim, lib_file, "advection"); - LibCoefficient mu(lib_file, "mu", false, mu_param); - LibCoefficient force(lib_file, "force"); - - // Define weak form and evolution - StabConvDifIntegrator integrator(adv, mu, force); - ParNonlinearForm form(space); - form.AddDomainIntegrator(&integrator); - form.UseExternalIntegrators(); - - Array ess_bdr(space->GetMesh()->bdr_attributes.Max()); - ess_bdr = 0; - for (int b = 0; b < strong_bdr.Size(); ++b) { - ess_bdr[strong_bdr[b]-1] = 1; + std::cout<<"Solve tau problem\n"; + // Define the gridfunction and solution vector + ParGridFunction phi_gf(space); + LibCoefficient sol_phi(lib_file, "sol_phi"); + phi_gf.ProjectCoefficient(sol_phi); + Vector xp; + phi_gf.GetTrueDofs(xp); + + // Define the visualisation output + VisItDataCollection vdc("tau", &pmesh); + vdc.SetPrefixPath(vis_dir); + vdc.RegisterField("phi", &phi_gf); + vdc.SetCycle(0); + vdc.Save(); + + // Define the physical parameters + LibVectorCoefficient adv(dim, lib_file, "advection"); + LibCoefficient mu(lib_file, "mu", false, mu_param); + LibCoefficient force(lib_file, "force"); + + // Define the inverse estimate + InverseEstimateCoefficient inv_est(space); + + // Define weak form and evolution + StabTauIntegrator integrator(adv, mu, inv_est); + ParNonlinearForm form(space); + form.AddDomainIntegrator(&integrator); + form.UseExternalIntegrators(); + + Array ess_bdr(space->GetMesh()->bdr_attributes.Max()); + ess_bdr = 0; + for (int b = 0; b < strong_bdr.Size(); ++b) + { + ess_bdr[strong_bdr[b]-1] = 1; + } + ess_bdr.Print(); + form.SetEssentialBC(ess_bdr); + //form.SetWeakBC (weak_bdr); + + Solver* pc_mom = nullptr; + HypreILU* ilu_mom = new HypreILU(); + pc_mom = ilu_mom; + + // Set up the Jacobian solver + FGMRESSolver gmres(MPI_COMM_WORLD); + gmres.iterative_mode = false; + gmres.SetRelTol(GMRES_RelTol); + gmres.SetMaxIter(GMRES_MaxIter); + gmres.SetKDim(GMRES_MaxIter+1); + gmres.SetPrintLevel(3); + gmres.SetPreconditioner(*pc_mom); + + // Set up the Newton solver + NewtonSolver newton_solver(MPI_COMM_WORLD); + newton_solver.SetOperator(form); + newton_solver.iterative_mode = true; + newton_solver.SetPrintLevel(1); + newton_solver.SetRelTol(Newton_RelTol); + newton_solver.SetMaxIter(Newton_MaxIter); + newton_solver.SetSolver(gmres); + + // Solver nonlinear system + Vector zero(space->TrueVSize()); + zero = 0.0; + newton_solver.Mult(zero, xp); + phi_gf.Distribute(xp); + + // Write solution + vdc.SetCycle(1); + vdc.Save(); } - ess_bdr.Print(); - form.SetEssentialBC(ess_bdr); - //form.SetWeakBC (weak_bdr); - - Solver* pc_mom = nullptr; - HypreILU* ilu_mom = new HypreILU(); - pc_mom = ilu_mom; - - // Set up the Jacobian solver - FGMRESSolver gmres(MPI_COMM_WORLD); - gmres.iterative_mode = false; - gmres.SetRelTol(GMRES_RelTol); - gmres.SetMaxIter(GMRES_MaxIter); - gmres.SetKDim(GMRES_MaxIter+1); - gmres.SetPrintLevel(3); - gmres.SetPreconditioner(*pc_mom); - - // Set up the Newton solver - NewtonSolver newton_solver(MPI_COMM_WORLD); - newton_solver.SetOperator(form); - newton_solver.iterative_mode = true; - newton_solver.SetPrintLevel(1); - newton_solver.SetRelTol(Newton_RelTol); - newton_solver.SetMaxIter(Newton_MaxIter); - newton_solver.SetSolver(gmres); - - // Solver nonlinear system - Vector zero(space->TrueVSize()); - zero = 0.0; - newton_solver.Mult(zero, xp); - phi_gf.Distribute(xp); - - // Compute errors - LibVectorCoefficient sol_grad(dim, lib_file, "grad_phi", false); - if (sol_grad.Foundfunction()) + + /* + + + CONDIF PART + + */ + { - int order_quad = max(2, 2*order+1); - const IntegrationRule *irs[Geometry::NumGeom]; - for (int i=0; i < Geometry::NumGeom; ++i) + std::cout<<"Solve convection-diffusion problem\n"; + // Define the gridfunction and solution vector + ParGridFunction phi_gf(space); + LibCoefficient sol_phi(lib_file, "sol_phi"); + phi_gf.ProjectCoefficient(sol_phi); + Vector xp; + phi_gf.GetTrueDofs(xp); + + // Define the visualisation output + VisItDataCollection vdc("phi", &pmesh); + vdc.SetPrefixPath(vis_dir); + vdc.RegisterField("phi", &phi_gf); + vdc.SetCycle(0); + vdc.Save(); + + // Define the physical parameters + LibVectorCoefficient adv(dim, lib_file, "advection"); + LibCoefficient mu(lib_file, "mu", false, mu_param); + LibCoefficient force(lib_file, "force"); + + // Define the inverse estimate + InverseEstimateCoefficient inv_est(space); + + // Define weak form and evolution + StabConvDifIntegrator integrator(adv, mu, force, inv_est); + ParNonlinearForm form(space); + form.AddDomainIntegrator(&integrator); + form.UseExternalIntegrators(); + + Array ess_bdr(space->GetMesh()->bdr_attributes.Max()); + ess_bdr = 0; + for (int b = 0; b < strong_bdr.Size(); ++b) { - irs[i] = &(IntRules.Get(i, order_quad)); + ess_bdr[strong_bdr[b]-1] = 1; } + ess_bdr.Print(); + form.SetEssentialBC(ess_bdr); + //form.SetWeakBC (weak_bdr); + + Solver* pc_mom = nullptr; + HypreILU* ilu_mom = new HypreILU(); + pc_mom = ilu_mom; + + // Set up the Jacobian solver + FGMRESSolver gmres(MPI_COMM_WORLD); + gmres.iterative_mode = false; + gmres.SetRelTol(GMRES_RelTol); + gmres.SetMaxIter(GMRES_MaxIter); + gmres.SetKDim(GMRES_MaxIter+1); + gmres.SetPrintLevel(3); + gmres.SetPreconditioner(*pc_mom); + + // Set up the Newton solver + NewtonSolver newton_solver(MPI_COMM_WORLD); + newton_solver.SetOperator(form); + newton_solver.iterative_mode = true; + newton_solver.SetPrintLevel(1); + newton_solver.SetRelTol(Newton_RelTol); + newton_solver.SetMaxIter(Newton_MaxIter); + newton_solver.SetSolver(gmres); + + // Solver nonlinear system + Vector zero(space->TrueVSize()); + zero = 0.0; + newton_solver.Mult(zero, xp); + phi_gf.Distribute(xp); + + // Compute errors + LibVectorCoefficient sol_grad(dim, lib_file, "grad_phi", false); + if (sol_grad.Foundfunction()) + { + int order_quad = max(2, 2*order+1); + const IntegrationRule *irs[Geometry::NumGeom]; + for (int i=0; i < Geometry::NumGeom; ++i) + { + irs[i] = &(IntRules.Get(i, order_quad)); + } - double err_phi = phi_gf.ComputeL2Error(sol_phi, irs); - double norm_phi = ComputeGlobalLpNorm(2., sol_phi, pmesh, irs); - std::cout << "|| phi_h - phi_ex || / || phi_ex || = " << err_phi / norm_phi << "\n"; + double err_phi = phi_gf.ComputeL2Error(sol_phi, irs); + double norm_phi = ComputeGlobalLpNorm(2., sol_phi, pmesh, irs); + std::cout << "|| phi_h - phi_ex || / || phi_ex || = " << err_phi / norm_phi << + "\n"; - err_phi = phi_gf.ComputeGradError(&sol_grad, irs); - norm_phi = ComputeGlobalLpNorm(2., sol_grad, pmesh, irs); - std::cout << "||grad phi_h - grad phi_ex || / || grad phi_ex || = " << err_phi / norm_phi << "\n"; - } + err_phi = phi_gf.ComputeGradError(&sol_grad, irs); + norm_phi = ComputeGlobalLpNorm(2., sol_grad, pmesh, irs); + std::cout << "||grad phi_h - grad phi_ex || / || grad phi_ex || = " << err_phi / + norm_phi << "\n"; + } - // Write solution - vdc.SetCycle(1); - vdc.Save(); + // Write solution + vdc.SetCycle(1); + vdc.Save(); - { - char vishost[] = "localhost"; - int visport = 19916; - socketstream sol_sock(vishost, visport); - sol_sock << "parallel " << num_procs << " " << myid << "\n"; - sol_sock.precision(8); - sol_sock << "solution\n" << pmesh << phi_gf << flush; + { + char vishost[] = "localhost"; + int visport = 19916; + socketstream sol_sock(vishost, visport); + sol_sock << "parallel " << num_procs << " " << myid << "\n"; + sol_sock.precision(8); + sol_sock << "solution\n" << pmesh << phi_gf << flush; + } } - // Free the used memory. delete fec; delete space; diff --git a/supg/tau_variance.cpp b/supg/tau_variance.cpp new file mode 100644 index 0000000..8fdd691 --- /dev/null +++ b/supg/tau_variance.cpp @@ -0,0 +1,259 @@ +// This file implements a small test to measure the variance of the stabilization +// parameter tau under permutations of element vertex ordering. The current +// method in supg/integrator.cpp uses GeomToPerfGeomJac, which should make tau +// invariant to vertex ordering for a single linear element because it maps to +// a perfect reference (equilateral triangle / regular tetrahedron, etc.). + +#include "mfem.hpp" +#include +#include +#include +#include +#include +#include +#include + +using namespace mfem; +using std::cout; +using std::endl; + +static std::string BuildSingleTriMeshText(const std::array,3> &v) +{ + // Reproduce mfem/data/ref-triangle.mesh layout but with custom vertex order + std::ostringstream os; + os << "MFEM mesh v1.0\n\n"; + os << "dimension\n2\n\n"; + os << "elements\n1\n"; + os << "1 2 0 1 2\n\n"; + os << "boundary\n3\n"; + os << "1 1 0 1\n"; + os << "2 1 1 2\n"; + os << "3 1 2 0\n\n"; + os << "vertices\n3\n2\n"; + os.setf(std::ios::fixed); os.precision(15); + for (int i = 0; i < 3; i++) + { + os << v[i][0] << ' ' << v[i][1] << "\n"; + } + os << "\n"; + return os.str(); +} + +static std::string BuildSingleTetMeshText(const std::array,4> &v) +{ + // Reproduce mfem/data/ref-tetrahedron.mesh layout but with custom vertex order + std::ostringstream os; + os << "MFEM mesh v1.0\n\n"; + os << "dimension\n3\n\n"; + os << "elements\n1\n"; + // Element: attribute=1, geometry=4 (TETRAHEDRON), vertices 0 1 2 3 + os << "1 4 0 1 2 3\n\n"; + // Boundary: 4 triangular faces with attributes 1..4 + os << "boundary\n4\n"; + os << "1 2 1 2 3\n"; + os << "2 2 0 3 2\n"; + os << "3 2 0 1 3\n"; + os << "4 2 0 2 1\n\n"; + // Vertices (4) with dimension flag 3 + os << "vertices\n4\n3\n"; + os.setf(std::ios::fixed); os.precision(15); + for (int i = 0; i < 4; i++) + { + os << v[i][0] << ' ' << v[i][1] << ' ' << v[i][2] << "\n"; + } + os << "\n"; + return os.str(); +} + +static double TauFromG(const DenseMatrix &G, const Vector &a, double k, double Ch) +{ + // Mirror supg/integrator.cpp: StabConvDifIntegrator::GetTau + // Production code does: Ch = 1/Ch; and then adds (Ch*Ch*k*k) inside the i,j loop + // i.e. acc = 1e-10 + sum_{i,j} ( G(i,j) a[i] a[j] + (1/Ch)^2 k^2 ) + double invCh = 1.0 / Ch; + double tau_acc = 1e-10; + const int dim = G.Width(); + for (int j = 0; j < dim; j++) + { + for (int i = 0; i < dim; i++) + { + tau_acc += G(i,j) * a[i] * a[j] + invCh * invCh * k * k; + } + } + return 1.0 / std::sqrt(tau_acc); +} + +static void ComputeTausForMeshText(const std::string &mesh_text, + double &tau_old, double &tau_new, + const Vector &a, double k, double Ch) +{ + std::istringstream is(mesh_text); + Mesh mesh(is, 1, 0, false); + + // Uncomment below to check for negative determinant + // int ninv = mesh.CheckElementOrientation(false); + // std::cout << "Inverted elements (without fixing): " << ninv << "\n"; + + // One element mesh, pick centroid of reference element (triangle/tetra) + ElementTransformation *T = mesh.GetElementTransformation(0); + IntegrationPoint ip; + if (mesh.Dimension() == 2) + { + ip.Set2(1.0/3.0, 1.0/3.0); + } + else if (mesh.Dimension() == 3) + { + ip.Set3(1.0/4.0, 1.0/4.0, 1.0/4.0); + } + else + { + MFEM_ABORT("Unsupported mesh dimension in tau_variance test."); + } + T->SetIntPoint(&ip); + + const DenseMatrix &invJ = T->InverseJacobian(); + + // Old method: G = invJ^T invJ + DenseMatrix G_old(invJ.Width()); + MultAtB(invJ, invJ, G_old); + tau_old = TauFromG(G_old, a, k, Ch); + + // New method using GeomToPerfGeomJac + const DenseMatrix &A = Geometries.GetGeomToPerfGeomJac(T->GetGeometryType()); + DenseMatrix invJ_perf(invJ.Height(), invJ.Width()); + Mult(A, invJ, invJ_perf); + DenseMatrix G_new(invJ.Width()); + MultAtB(invJ_perf, invJ_perf, G_new); + tau_new = TauFromG(G_new, a, k, Ch); +} + +// -------------------- Small utilities and geometry-specific runners -------------------- + +struct Stats +{ + double mean = 0.0; + double var = 0.0; + double mn = 0.0; + double mx = 0.0; +}; + +static Stats ComputeStats(const std::vector &x) +{ + MFEM_VERIFY(!x.empty(), "ComputeStats requires non-empty input"); + Stats s{}; + for (double v : x) s.mean += v; + s.mean /= x.size(); + for (double v : x) { const double d = v - s.mean; s.var += d*d; } + s.var /= x.size(); + s.mn = x[0]; s.mx = x[0]; + for (double v : x) { if (v < s.mn) s.mn = v; if (v > s.mx) s.mx = v; } + return s; +} + +template +static void PrintPermutationValues(const std::vector> &perms, + const std::vector &taus_old, + const std::vector &taus_new) +{ + using std::cout; using std::endl; + cout << "\nValues per permutation ("; + for (size_t i = 0; i < N; i++) { cout << 'p' << i << (i+1 +#include +#include +#include + +using namespace mfem; +using std::cout; +using std::endl; + + +/// Create a mesh with one triangle, with different element vertices +static Mesh *BuildSingleTriMesh(const int shift) +{ + + // Reproduce mfem/data/ref-triangle.mesh layout but with custom vertex order + std::ostringstream os; + os << "MFEM mesh v1.0\n\n"; + os << "dimension\n2\n\n"; + os << "elements\n1\n"; + os << "1 2"; + for (int i = 0; i < 3; i++) + { + os << " "<<(i + shift)%3; + } + os << "\n\n"; + os << "boundary\n3\n"; + os << "1 1 0 1\n"; + os << "2 1 1 2\n"; + os << "3 1 2 0\n\n"; + os << "vertices\n3\n2\n"; + os << "0 0\n"; + os << "1 0\n"; + os << "0 1\n"; + os << "\n"; + std::istringstream is(os.str()); + return new Mesh(is, 1, 0, false); +} + +/// Create a mesh with one tetrahadron, with different element vertices +static Mesh *BuildSingleTetMesh(const int shift) +{ + // Reproduce mfem/data/ref-tetrahedron.mesh layout but with custom vertex order + std::ostringstream os; + os << "MFEM mesh v1.0\n\n"; + os << "dimension\n3\n\n"; + os << "elements\n1\n"; + // Element: attribute=1, geometry=4 (TETRAHEDRON), vertices 0 1 2 3 + os << "1 4"; + for (int i = 0; i < 4; i++) + { + os << " "<<(i + shift)%4; + } + os << "\n\n"; + // Boundary: 4 triangular faces with attributes 1..4 + os << "boundary\n4\n"; + os << "1 2 1 2 3\n"; + os << "2 2 0 3 2\n"; + os << "3 2 0 1 3\n"; + os << "4 2 0 2 1\n\n"; + // Vertices (4) with dimension flag 3 + os << "vertices\n4\n3\n"; + os << "0 0 0\n"; + os << "1 0 0\n"; + os << "0 1 0\n"; + os << "0 0 1\n"; + os << "\n"; + std::istringstream is(os.str()); + return new Mesh(is, 1, 0, false); +} + +static double TauFromG(const DenseMatrix &G, const Vector &a) +{ + double tau_acc = 1e-10; + const int dim = G.Width(); + for (int j = 0; j < dim; j++) + { + for (int i = 0; i < dim; i++) + { + tau_acc += G(i,j) * a[i] * a[j]; + } + } + return 1.0 / std::sqrt(tau_acc); +} + +static void ComputeConvectiveTau(Mesh *mesh, + double &tau_old, double &tau_new, + const Vector &a) +{ + // Uncomment below to check for negative determinant + // int ninv = mesh.CheckElementOrientation(false); + // std::cout << "Inverted elements (without fixing): " << ninv << "\n"; + + // One element mesh, pick a canonical interior point for the reference element + ElementTransformation *T = mesh->GetElementTransformation(0); + IntegrationPoint ip; + switch (T->GetGeometryType()) + { + case Geometry::TRIANGLE: ip.Set2(1.0/3.0, 1.0/3.0); break; // centroid + case Geometry::TETRAHEDRON: ip.Set3(1.0/4.0, 1.0/4.0, 1.0/4.0); break; // centroid + case Geometry::PRISM: ip.Set3(1.0/3.0, 1.0/3.0, 0.5); break; // tri-centroid x mid-height + case Geometry::PYRAMID: ip.Set3(0.375, 0.375, 0.25); break; // ref pyramid centroid? + default: MFEM_ABORT("Unsupported geometry in tau_variance_alt test."); + } + T->SetIntPoint(&ip); + + const DenseMatrix &invJ = T->InverseJacobian(); + + // Old method: G = invJ^T invJ + DenseMatrix G_old(invJ.Width()); + MultAtB(invJ, invJ, G_old); + tau_old = TauFromG(G_old, a); + + // New method using GeomToPerfGeomJac + const DenseMatrix &A = Geometries.GetGeomToPerfGeomJac(T->GetGeometryType()); + DenseMatrix invJ_perf(invJ.Height(), invJ.Width()); + Mult(A, invJ, invJ_perf); + DenseMatrix G_new(invJ.Width()); + MultAtB(invJ_perf, invJ_perf, G_new); + tau_new = TauFromG(G_new, a); +} + +static void ComputeDiffusive(Mesh *mesh, + double &tr_old, double &tr_new, + double &nf_old, double &nf_new) +{ + // One element mesh, pick a canonical interior point for the reference element + ElementTransformation *T = mesh->GetElementTransformation(0); + IntegrationPoint ip; + switch (T->GetGeometryType()) + { + case Geometry::TRIANGLE: ip.Set2(1.0/3.0, 1.0/3.0); break; + case Geometry::TETRAHEDRON: ip.Set3(1.0/4.0, 1.0/4.0, 1.0/4.0); break; + case Geometry::PRISM: ip.Set3(1.0/3.0, 1.0/3.0, 0.5); break; + case Geometry::PYRAMID: ip.Set3(0.375, 0.375, 0.25); break; + default: MFEM_ABORT("Unsupported geometry in tau_variance_alt test."); + } + T->SetIntPoint(&ip); + + const DenseMatrix &invJ = T->InverseJacobian(); + + // Old method: G = invJ^T invJ + DenseMatrix G_old(invJ.Width()); + MultAtB(invJ, invJ, G_old); + tr_old = G_old.Trace(); + nf_old = G_old.FNorm(); + + // New method using GeomToPerfGeomJac + const DenseMatrix &A = Geometries.GetGeomToPerfGeomJac(T->GetGeometryType()); + DenseMatrix invJ_perf(invJ.Height(), invJ.Width()); + Mult(A, invJ, invJ_perf); + DenseMatrix G_new(invJ.Width()); + MultAtB(invJ_perf, invJ_perf, G_new); + tr_new = G_new.Trace(); + nf_new = G_new.FNorm(); +} + + +// -------------------- Small utilities and geometry-specific runners -------------------- + +struct Stats +{ + double mean = 0.0; + double var = 0.0; + double mn = 0.0; + double mx = 0.0; +}; + +static Stats ComputeStats(const std::vector &x) +{ + MFEM_VERIFY(!x.empty(), "ComputeStats requires non-empty input"); + Stats s{}; + for (double v : x) s.mean += v; + s.mean /= x.size(); + for (double v : x) { const double d = v - s.mean; s.var += d*d; } + s.var /= x.size(); + s.mn = x[0]; s.mx = x[0]; + for (double v : x) { if (v < s.mn) s.mn = v; if (v > s.mx) s.mx = v; } + return s; +} + + +static void PrintPermutationValues(const std::vector &taus_old, + const std::vector &taus_new) +{ + for (size_t i = 0; i < taus_old.size(); i++) + { + cout << "perm ("; + for (size_t j = 0; j < taus_old.size(); j++) + { + cout <<" "<< (i+j)%taus_old.size(); + } + cout << " )"; + cout << " : old=" << taus_old[i] << ", new=" << taus_new[i] << endl; + } +} + +static int RunTriangleTest() +{ + Mesh *mesh[3]; + for (size_t i = 0; i < 3; i++) + { + mesh[i] = BuildSingleTriMesh(i); + } + + int num_samples = 101; + Vector a(2); + std::ofstream out ("conv_tau_2D.dat"); + for (size_t i = 0; i < num_samples; i++) + { + double theta = (2*M_PI*i)/(num_samples-1); + a[0] = sin(theta); + a[1] = cos(theta); + out<< theta ; + double to, tn; + for (size_t i = 0; i < 3; i++) + { + ComputeConvectiveTau(mesh[i], to, tn, a); + out<< " "< trs_old, trs_new; + std::vector nfs_old, nfs_new; + double to, tn, no, nn; + for (size_t i = 0; i < 3; i++) + { + ComputeDiffusive(mesh[i], to, tn, no, nn); + trs_old.push_back(to); + trs_new.push_back(tn); + nfs_old.push_back(no); + nfs_new.push_back(nn); + } + + Stats sto = ComputeStats(trs_old); + Stats stn = ComputeStats(trs_new); + Stats sno = ComputeStats(nfs_old); + Stats snn = ComputeStats(nfs_new); + + cout << "Trace variance test over all 3 permutations of a single triangle" << endl; + + cout.setf(std::ios::scientific); cout.precision(6); + cout << "Old method (raw invJ): mean=" << sto.mean << ", var=" << sto.var + << ", min=" << sto.mn << ", max=" << sto.mx << endl; + cout << "New method (perf geom): mean=" << stn.mean << ", var=" << stn.var + << ", min=" << stn.mn << ", max=" << stn.mx << endl; + + PrintPermutationValues(trs_old, trs_new); + + cout << "FNorm variance test over all 3 permutations of a single triangle" << endl; + cout.setf(std::ios::scientific); cout.precision(6); + cout << "Old method (raw invJ): mean=" << sno.mean << ", var=" << sno.var + << ", min=" << sno.mn << ", max=" << sno.mx << endl; + cout << "New method (perf geom): mean=" << snn.mean << ", var=" << snn.var + << ", min=" << snn.mn << ", max=" << snn.mx << endl; + + PrintPermutationValues(nfs_old, nfs_new); + + for (size_t i = 0; i < 3; i++) + { + delete mesh[i]; + } + // Same failure criterion used previously + return 0;//(sn.var > 1e-28 && sn.var > 1e-6 * so.var) ? 1 : 0; +} + +static int RunTetrahedronTest() +{ + Mesh *mesh[4]; + for (size_t i = 0; i < 4; i++) + { + mesh[i] = BuildSingleTetMesh(i); + } + + int num_samples = 101; + Vector a(3); + + std::ofstream out ("conv_tau_3D.dat"); + for (size_t ii = 0; ii < num_samples; ii++) + { + double phi = (2*M_PI*ii)/(num_samples-1); + for (size_t i = 0; i < num_samples; i++) + { + double theta = (2*M_PI*i)/(num_samples-1); + + a[0] = sin(theta); + a[1] = cos(theta)*sin(phi); + a[2] = cos(theta)*cos(phi); + out<< theta <<" " << phi; + double to, tn; + for (size_t i = 0; i < 4; i++) + { + ComputeConvectiveTau(mesh[i], to, tn, a); + out<< " "< trs_old, trs_new; + std::vector nfs_old, nfs_new; + + double to, tn, no, nn; + for (size_t i = 0; i < 4; i++) + { + ComputeDiffusive(mesh[i], to, tn, no, nn); + trs_old.push_back(to); + trs_new.push_back(tn); + nfs_old.push_back(no); + nfs_new.push_back(nn); + } + + Stats sto = ComputeStats(trs_old); + Stats stn = ComputeStats(trs_new); + Stats sno = ComputeStats(nfs_old); + Stats snn = ComputeStats(nfs_new); + + + cout << "Trace variance test over all 4 permutations of a single tetrahedron" << endl; + + cout.setf(std::ios::scientific); cout.precision(6); + cout << "Old method (raw invJ): mean=" << sto.mean << ", var=" << sto.var + << ", min=" << sto.mn << ", max=" << sto.mx << endl; + cout << "New method (perf geom): mean=" << stn.mean << ", var=" << stn.var + << ", min=" << stn.mn << ", max=" << stn.mx << endl; + + PrintPermutationValues(trs_old, trs_new); + + cout << "FNorm variance test over all 4 permutations of a single tetrahedron" << endl; + cout.setf(std::ios::scientific); cout.precision(6); + cout << "Old method (raw invJ): mean=" << sno.mean << ", var=" << sno.var + << ", min=" << sno.mn << ", max=" << sno.mx << endl; + cout << "New method (perf geom): mean=" << snn.mean << ", var=" << snn.var + << ", min=" << snn.mn << ", max=" << snn.mx << endl; + + PrintPermutationValues(nfs_old, nfs_new); + + for (size_t i = 0; i < 4; i++) + { + delete mesh[i]; + } + + return 0; +} + +/// Create a mesh with one right prism, with different element vertex permutations +static Mesh *BuildSinglePrismMesh(const int rot) +{ + std::ostringstream os; + os << "MFEM mesh v1.0\n\n"; + os << "dimension\n3\n\n"; + os << "elements\n1\n"; + // Prism local vertex order: bottom tri (0,1,2), top tri (3,4,5) + // We only rotate cyclically the (0,1,2) and (3,4,5) blocks by 'rot' in [0,2] + int bot[3] = {0,1,2}; + int top[3] = {3,4,5}; + int b0 = bot[(0+rot)%3], b1 = bot[(1+rot)%3], b2 = bot[(2+rot)%3]; + int t0 = top[(0+rot)%3], t1 = top[(1+rot)%3], t2 = top[(2+rot)%3]; + os << "1 6 " << b0 << " " << b1 << " " << b2 << " " << t0 << " " << t1 << " " << t2 << "\n\n"; + os << "boundary\n5\n"; + // bottom triangle + os << "1 2 " << b0 << " " << b2 << " " << b1 << "\n"; + // top triangle + os << "1 2 " << t0 << " " << t1 << " " << t2 << "\n"; + // side quads + os << "1 3 " << b0 << " " << b1 << " " << t1 << " " << t0 << "\n"; + os << "1 3 " << b1 << " " << b2 << " " << t2 << " " << t1 << "\n"; + os << "1 3 " << b2 << " " << b0 << " " << t0 << " " << t2 << "\n\n"; + os << "vertices\n6\n3\n"; + os << "0 0 0\n"; + os << "1 0 0\n"; + os << "0 1 0\n"; + os << "0 0 1\n"; + os << "1 0 1\n"; + os << "0 1 1\n\n"; + std::istringstream is(os.str()); + return new Mesh(is, 1, 0, false); +} + +/// Create a mesh with one right pyramid, with different element vertex permutations +static Mesh *BuildSinglePyramidMesh(const int rot) +{ + std::ostringstream os; + os << "MFEM mesh v1.0\n\n"; + os << "dimension\n3\n\n"; + os << "elements\n1\n"; + // Pyramid local order: square base (0,1,2,3) + apex (4) + // Permute by cyclic rotation of base only by 'rot' in [0,3] + int base[4] = {0,1,2,3}; + int b0 = base[(0+rot)%4], b1 = base[(1+rot)%4], b2 = base[(2+rot)%4], b3 = base[(3+rot)%4]; + int apex = 4; + os << "1 7 " << b0 << " " << b1 << " " << b2 << " " << b3 << " " << apex << "\n\n"; + os << "boundary\n5\n"; + // base quad + os << "1 3 " << b3 << " " << b2 << " " << b1 << " " << b0 << "\n"; + // side triangles + os << "1 2 " << b0 << " " << b1 << " " << apex << "\n"; + os << "1 2 " << b1 << " " << b2 << " " << apex << "\n"; + os << "1 2 " << b2 << " " << b3 << " " << apex << "\n"; + os << "1 2 " << b3 << " " << b0 << " " << apex << "\n\n"; + os << "vertices\n5\n3\n"; + os << "0 0 0\n"; + os << "1 0 0\n"; + os << "1 1 0\n"; + os << "0 1 0\n"; + os << "0 0 1\n\n"; + std::istringstream is(os.str()); + return new Mesh(is, 1, 0, false); +} + +static int RunPrismTest() +{ + Mesh *mesh[3]; + for (size_t i = 0; i < 3; i++) { mesh[i] = BuildSinglePrismMesh((int)i); } + + int num_samples = 101; + Vector a(3); + + std::ofstream out("conv_tau_prism.dat"); + for (size_t ii = 0; ii < num_samples; ii++) + { + double phi = (2*M_PI*ii)/(num_samples-1); + for (size_t i = 0; i < num_samples; i++) + { + double theta = (2*M_PI*i)/(num_samples-1); + a[0] = sin(theta); + a[1] = cos(theta)*sin(phi); + a[2] = cos(theta)*cos(phi); + out << theta << " " << phi; + double to, tn; + for (size_t p = 0; p < 3; p++) + { + ComputeConvectiveTau(mesh[p], to, tn, a); + out << " " << to << " " << tn; + } + out << "\n"; + } + out << "\n"; + } + + std::vector trs_old, trs_new; + std::vector nfs_old, nfs_new; + double to, tn, no, nn; + for (size_t i = 0; i < 3; i++) + { + ComputeDiffusive(mesh[i], to, tn, no, nn); + trs_old.push_back(to); trs_new.push_back(tn); + nfs_old.push_back(no); nfs_new.push_back(nn); + } + + Stats sto = ComputeStats(trs_old); + Stats stn = ComputeStats(trs_new); + Stats sno = ComputeStats(nfs_old); + Stats snn = ComputeStats(nfs_new); + + cout << "Trace variance test over 3 cyclic base-rotations of a single prism" << endl; + cout.setf(std::ios::scientific); cout.precision(6); + cout << "Old method (raw invJ): mean=" << sto.mean << ", var=" << sto.var + << ", min=" << sto.mn << ", max=" << sto.mx << endl; + cout << "New method (perf geom): mean=" << stn.mean << ", var=" << stn.var + << ", min=" << stn.mn << ", max=" << stn.mx << endl; + PrintPermutationValues(trs_old, trs_new); + + cout << "FNorm variance test over 3 cyclic base-rotations of a single prism" << endl; + cout.setf(std::ios::scientific); cout.precision(6); + cout << "Old method (raw invJ): mean=" << sno.mean << ", var=" << sno.var + << ", min=" << sno.mn << ", max=" << sno.mx << endl; + cout << "New method (perf geom): mean=" << snn.mean << ", var=" << snn.var + << ", min=" << snn.mn << ", max=" << snn.mx << endl; + PrintPermutationValues(nfs_old, nfs_new); + + for (size_t i = 0; i < 3; i++) { delete mesh[i]; } + return 0; +} + +static int RunPyramidTest() +{ + Mesh *mesh[4]; + for (size_t i = 0; i < 4; i++) { mesh[i] = BuildSinglePyramidMesh((int)i); } + + int num_samples = 101; + Vector a(3); + + std::ofstream out("conv_tau_pyramid.dat"); + for (size_t ii = 0; ii < num_samples; ii++) + { + double phi = (2*M_PI*ii)/(num_samples-1); + for (size_t i = 0; i < num_samples; i++) + { + double theta = (2*M_PI*i)/(num_samples-1); + a[0] = sin(theta); + a[1] = cos(theta)*sin(phi); + a[2] = cos(theta)*cos(phi); + out << theta << " " << phi; + double to, tn; + for (size_t p = 0; p < 4; p++) + { + ComputeConvectiveTau(mesh[p], to, tn, a); + out << " " << to << " " << tn; + } + out << "\n"; + } + out << "\n"; + } + + std::vector trs_old, trs_new; + std::vector nfs_old, nfs_new; + double to, tn, no, nn; + for (size_t i = 0; i < 4; i++) + { + ComputeDiffusive(mesh[i], to, tn, no, nn); + trs_old.push_back(to); trs_new.push_back(tn); + nfs_old.push_back(no); nfs_new.push_back(nn); + } + + Stats sto = ComputeStats(trs_old); + Stats stn = ComputeStats(trs_new); + Stats sno = ComputeStats(nfs_old); + Stats snn = ComputeStats(nfs_new); + + cout << "Trace variance test over 4 cyclic base-rotations of a single pyramid" << endl; + cout.setf(std::ios::scientific); cout.precision(6); + cout << "Old method (raw invJ): mean=" << sto.mean << ", var=" << sto.var + << ", min=" << sto.mn << ", max=" << sto.mx << endl; + cout << "New method (perf geom): mean=" << stn.mean << ", var=" << stn.var + << ", min=" << stn.mn << ", max=" << stn.mx << endl; + PrintPermutationValues(trs_old, trs_new); + + cout << "FNorm variance test over 4 cyclic base-rotations of a single pyramid" << endl; + cout.setf(std::ios::scientific); cout.precision(6); + cout << "Old method (raw invJ): mean=" << sno.mean << ", var=" << sno.var + << ", min=" << sno.mn << ", max=" << sno.mx << endl; + cout << "New method (perf geom): mean=" << snn.mean << ", var=" << snn.var + << ", min=" << snn.mn << ", max=" << snn.mx << endl; + PrintPermutationValues(nfs_old, nfs_new); + + for (size_t i = 0; i < 4; i++) { delete mesh[i]; } + return 0; +} + +int main(int argc, char** argv) +{ + int fail = RunTriangleTest(); + fail |= RunTetrahedronTest(); + fail |= RunPrismTest(); + fail |= RunPyramidTest(); + + return fail; +} diff --git a/util/coefficients.cpp b/util/coefficients.cpp index cec9e2b..e6e2ee9 100644 --- a/util/coefficients.cpp +++ b/util/coefficients.cpp @@ -11,6 +11,206 @@ using namespace std; using namespace mfem; +extern "C" void +dsygvx_(int *ITYPE, char *JOBZ, char *RANGE, char *UPLO, + int *N, double *A, int *LDA, double *B, int *LDB, + double *VL, double *VU, int *IL, int *IU, + double *ABSTOL, int *M, double *W, double *Z, + int *LDZ, double *WORK, int *LWORK,int *IWORK, + int *IFAIL, int *INFO); + +void dsygvx_Eigensystem(DenseMatrix &a, DenseMatrix &b, + Vector &ev, DenseMatrix *evect, + char RANGE, real_t VL, real_t VU, int IL, int IU) +{ +#ifdef MFEM_USE_LAPACK + ev.SetSize(a.Width()); + int ITYPE = 1; + char JOBZ = 'N'; + char UPLO = 'U'; + int N = a.Width(); + real_t *A = new real_t[N*N]; + int LDA = N; + real_t *B = new real_t[N*N]; + int LDB = N; + real_t ABSTOL = 0.0; + int M; + real_t *W = ev.GetData(); + real_t *Z = NULL; + int LDZ = 1; + real_t QWORK; + real_t *WORK = NULL; + int LWORK = -1; // query optimal (double) workspace size + int *IWORK = new int[5*N]; + int *IFAIL = new int[N]; + int INFO; + + if (evect) // Compute eigenvectors too + { + evect->SetSize(N); + + JOBZ = 'V'; + Z = evect->Data(); + LDZ = N; + } + + int hw = a.Height() * a.Width(); + real_t *data_a = a.Data(); + real_t *data_b = b.Data(); + + for (int i = 0; i < hw; i++) + { + A[i] = data_a[i]; + B[i] = data_b[i]; + } + +// MFEM_LAPACK_PREFIX() + dsygvx_( &ITYPE, &JOBZ, &RANGE, &UPLO, + &N, A, &LDA, B, &LDB, &VL, &VU, &IL, &IU, + &ABSTOL, &M, W, Z, &LDZ, &QWORK, &LWORK, + IWORK, IFAIL, &INFO ); + + LWORK = (int) QWORK; + WORK = new real_t[LWORK]; + + //MFEM_LAPACK_PREFIX() + dsygvx_( &ITYPE, &JOBZ, &RANGE, &UPLO, + &N, A, &LDA, B, &LDB, &VL, &VU, &IL, &IU, + &ABSTOL, &M, W, Z, &LDZ, WORK, &LWORK, + IWORK, IFAIL, &INFO ); + + if (INFO != 0) + { + mfem::err << "dsyevr_Eigensystem(...): DSYEVR error code: " + << INFO << endl; + mfem_error(); + } + + if (evect) // Compute eigenvectors too + { + evect->SetSize(N,M); + } + delete [] IFAIL; + delete [] IWORK; + delete [] WORK; + delete [] B; + if (evect == NULL) { delete [] A; } + +#else + MFEM_CONTRACT_VAR(a); + MFEM_CONTRACT_VAR(ev); + MFEM_CONTRACT_VAR(evect); +#endif +} + +real_t Eigenvalue(DenseMatrix &a, DenseMatrix &b, int i = -1) +{ +#ifdef MFEM_USE_LAPACK + if (i < 0) { i = a.Width(); } + Vector ev; + dsygvx_Eigensystem(a, b, ev, NULL, 'I', 0.0, 0.0, i, i); + return ev[0]; + +#else + MFEM_CONTRACT_VAR(ns); + MFEM_CONTRACT_VAR(tol); + mfem_error("DenseMatrix::Eigenvalue: Compiled without LAPACK"); + return 0.0; +#endif +} + +InverseEstimateCoefficient::InverseEstimateCoefficient(FiniteElementSpace *f) + : fes(f), Q(NULL), ir(NULL) +{ + ComputeInverseEstimates(); +} + +InverseEstimateCoefficient::InverseEstimateCoefficient(FiniteElementSpace *f, + Coefficient &q) + : fes(f), Q(&q), ir(NULL) +{ + ComputeInverseEstimates(); +} + +GridFunction *InverseEstimateCoefficient::GetGridFunction() +{ + FiniteElementCollection* fec_ec = new L2_FECollection(0, + fes ->GetMesh()->Dimension()); + FiniteElementSpace *fes_ec = new FiniteElementSpace(fes ->GetMesh(), fec_ec); + GridFunction *gf = new GridFunction(fes_ec, elemInvEst.GetData()); + gf->MakeOwner(fec_ec); + return gf; +} + +void InverseEstimateCoefficient::ComputeInverseEstimates() +{ + elemInvEst.SetSize(fes -> GetNE()); + SetIntRule(*fes->GetFE(0)); + for (int i = 0; i < fes -> GetNE(); i++) + { + elemInvEst[i] = ElementInverseEstimate(*fes->GetFE(i), + *fes->GetElementTransformation(i)); + } +} + +void InverseEstimateCoefficient::SetIntRule(const FiniteElement &el) +{ + ir = &IntRules.Get(el.GetGeomType(), 2*el.GetOrder()); +} + +real_t InverseEstimateCoefficient::ElementInverseEstimate( + const FiniteElement &el, + ElementTransformation &Trans) +{ + if (el.GetOrder() < 2) + { + return std::numeric_limits::min(); + } + + int nd = el.GetDof(); + int dim = el.GetDim(); + + shape.SetSize(nd); + dshape.SetSize(nd,dim); + laplace.SetSize(nd); + + lapmat.SetSize(nd,nd); + bimat.SetSize(nd,nd); + ovec.SetSize(nd); + + real_t w,q = 1.0; + + bimat = 0.0; + lapmat = 0.0; + ovec = 0.0; + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + Trans.SetIntPoint(&ip); + w = Trans.Weight()*ip.weight; + if (Q) + { + q = Q->Eval(Trans, ip); + } + + el.CalcPhysDShape(Trans, dshape); + AddMult_a_AAt(w*q, dshape, lapmat); + + el.CalcPhysLaplacian(Trans, laplace); + AddMult_a_VVt(w*q*q, laplace, bimat); + + el.CalcPhysShape(Trans, shape); + ovec.Add(w, shape); + } + ovec *= 1.0/ovec.Norml2(); + + // Correct nullspace + AddMultVVt(ovec, lapmat); + + // Return largest eigenvalue + return Eigenvalue(bimat, lapmat); +} + // Get the function from the library void LibCoefficient::GetLibFunction(string libName, vector funNames, diff --git a/util/coefficients.hpp b/util/coefficients.hpp index 341be32..e3b4c75 100644 --- a/util/coefficients.hpp +++ b/util/coefficients.hpp @@ -13,6 +13,61 @@ using namespace mfem; using namespace std; + + +class InverseEstimateCoefficient : public Coefficient +{ +private: + /// + Vector elemInvEst; + /// FE space on which the grid function lives. Owned if #fec is not NULL. + FiniteElementSpace *fes; + + /// + const IntegrationRule *ir; + + /// + Coefficient *Q; + Vector laplace, shape, ovec, evec; + DenseMatrix dshape, lapmat, bimat; + + /// + void SetIntRule(const FiniteElement &el); + + /// + void ComputeInverseEstimates(); + + real_t ElementInverseEstimate(const FiniteElement &el, + ElementTransformation &Trans); + +public: + /// + InverseEstimateCoefficient(FiniteElementSpace *f); + InverseEstimateCoefficient(FiniteElementSpace *f, Coefficient &q); + + /// Caller gets owner ship of GridFunction and + GridFunction *GetGridFunction(); + + /// Reset the scalar factor + void SetDiffusion(Coefficient &q) + { + if (Q != &q) + { + Q = &q; + ComputeInverseEstimates(); + } + } + /// Return the scalar factor + Coefficient * GetDiffusion() const { return Q; } + + /// Evaluate the coefficient at @a ip. + virtual real_t Eval(ElementTransformation &T, + const IntegrationPoint &ip) + { return elemInvEst[T.ElementNo]; } + +}; + + /// Coefficient class with a function defined in a seperate C-function class LibCoefficient : public Coefficient { @@ -63,7 +118,7 @@ class LibCoefficient : public Coefficient bool Foundfunction(){ return TDFunction;} - /// Evaluate + /// Evaluate using Coefficient::Eval; real_t Eval(ElementTransformation &T,