FreeSurfaceOps.hpp

Go to the documentation of this file.

/**
 * @file FreeSurfaceOps.hpp
 * @brief Implements operators for assembling residuals and Jacobians for the Navier–Stokes–Cahn–Hilliard free-surface problem.
 * 
 * This file contains boundary and domain operators used in MoFEM pipelines:
 * - Boundary operators: normal constraints, wetting angle, lift force.
 * - Domain operators: momentum, phase evolution, chemical potential.
 * 
 * Documentation: see `free_surface.dox` in the `doc` folder of this tutorial.
 * Note that documentation has been produced with the help of copilot. (Though all has had human review and editing!)
 * 
 * @see free_surface.dox
 * @see Lovrić et al., "Low Order Finite Element Methods for the Navier–Stokes–Cahn–Hilliard Equations", arXiv:1911.06718
 */
 
namespace FreeSurfaceOps {
 
 
/**
 * @brief Calculate lift force on free surface boundary
 * @param field_name Name of the field associated with the operator
 * @param p_ptr Pointer to the pressure vector on the free surface
 * @param lift_ptr Pointer to the lift force vector to be computed
 * @param ents_ptr Pointer to the range of entities (edges) to consider
 * @return MoFEMErrorCode
 * 
 * Computes the integral ∫_Γ ​-p n dΓ, where:
 * - p : pressure on the free surface
 * - n : outward normal vector on the boundary
 */
struct OpCalculateLift : public BoundaryEleOp {
  OpCalculateLift(const std::string field_name,
                  boost::shared_ptr<VectorDouble> p_ptr,
                  boost::shared_ptr<VectorDouble> lift_ptr,
                  boost::shared_ptr<Range> ents_ptr)
      : BoundaryEleOp(field_name, field_name, BoundaryEleOp::OPROW),
        pPtr(p_ptr), liftPtr(lift_ptr), entsPtr(ents_ptr) {
    std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
    doEntities[MBEDGE] = true;
  }
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    const auto fe_ent = getFEEntityHandle();
    if (entsPtr->find(fe_ent) != entsPtr->end()) {
 
      auto t_w = getFTensor0IntegrationWeight();
      auto t_p = getFTensor0FromVec(*pPtr);
      auto t_normal = getFTensor1Normal();
      auto t_coords = getFTensor1CoordsAtGaussPts();
      auto t_lift = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(*liftPtr);
 
      const auto nb_int_points = getGaussPts().size2();
 
      for (int gg = 0; gg != nb_int_points; gg++) {
 
        const double r = t_coords(0);
        const double alpha = cylindrical(r) * t_w;
        t_lift(i) -= t_normal(i) * (t_p * alpha); // −p n
 
        ++t_w;
        ++t_p;
        ++t_coords;
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<VectorDouble> pPtr;
  boost::shared_ptr<VectorDouble> liftPtr;
  boost::shared_ptr<Range> entsPtr;
};
 
//! [OpNormalConstrainRhs]
/**
 * @brief Boundary normal constraint RHS
 * @param field_name Name of the field associated with the operator
 * @param u_ptr Pointer to the velocity matrix on the free surface
 * @return MoFEMErrorCode
 * 
 * Computes the integral ∫_Γ ​w_L​·(n·u) dΓ, where:
 * - w_L : test function for Lagrange multiplier (λ)
 * - n   : outward normal vector on the boundary
 */
struct OpNormalConstrainRhs : public AssemblyBoundaryEleOp {
  OpNormalConstrainRhs(const std::string field_name,
                       boost::shared_ptr<MatrixDouble> u_ptr)
      : AssemblyBoundaryEleOp(field_name, field_name,
                              AssemblyBoundaryEleOp::OPROW),
        uPtr(u_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
    MoFEMFunctionBegin;
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_normal = getFTensor1Normal();
    auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
    auto t_row_base = row_data.getFTensor0N();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
      const double r = t_coords(0);
      const double alpha = t_w * cylindrical(r);
 
      int bb = 0;
      for (; bb != nbRows; ++bb) {
        locF[bb] += alpha * t_row_base * (t_normal(i) * t_u(i)); // λ * (n · u)
        ++t_row_base;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb)
        ++t_row_base;
 
      ++t_w;
      ++t_u;
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> uPtr;
};
//! [OpNormalConstrainRhs]
 
/**
 * @brief Boundary normal force RHS
 * @param field_name Name of the field associated with the operator
 * @param lambda_ptr Pointer to the Lagrange multiplier vector on the free surface
 * @return MoFEMErrorCode
 * 
 * Assembles the contribution of the Lagrange multiplier (λ) normal force
 * ie ∫_Γ N_u · (λ n) dΓ into the momentum residual.
 */
struct OpNormalForceRhs : public AssemblyBoundaryEleOp {
  OpNormalForceRhs(const std::string field_name,
                   boost::shared_ptr<VectorDouble> lambda_ptr)
      : AssemblyBoundaryEleOp(field_name, field_name,
                              AssemblyDomainEleOp::OPROW),
        lambdaPtr(lambda_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
    MoFEMFunctionBegin;
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_normal = getFTensor1Normal();
    auto t_lambda = getFTensor0FromVec(*lambdaPtr);
    auto t_row_base = row_data.getFTensor0N();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
      auto t_nf = getFTensor1FromArray<U_FIELD_DIM, U_FIELD_DIM>(locF);
 
      const double r = t_coords(0);
      const double alpha = t_w * cylindrical(r);
 
      int bb = 0;
      for (; bb != nbRows / U_FIELD_DIM; ++bb) {
 
        t_nf(i) += alpha * t_row_base * t_normal(i) * t_lambda; // N_u · (λ n)
        ++t_row_base;
        ++t_nf;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb)
        ++t_row_base;
 
      ++t_w;
      ++t_lambda;
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<VectorDouble> lambdaPtr;
};
 
//! [OpWettingAngleRhs]
/**
 * @brief Boundary wetting-angle RHS
 * @param field_name Name of the field associated with the operator
 * @param grad_h_ptr Pointer to the gradient of the free surface height matrix
 * @param ents_ptr Pointer to the range of entities (edges) to consider
 * @param wetting_angle Contact wetting angle in degrees
 * @return MoFEMErrorCode
 * 
 * Implements the boundary term from eq. 3.2d:
 * - ∫_Γ sh ε² ||∇ϕ|| cos(α) dΓ
 * The operator computes rhs_wetting = s * η² * ||∇h|| * cos(angle) and
 * assembles it into the G residual on selected meshset edges.
 */
struct OpWettingAngleRhs : public AssemblyBoundaryEleOp {
  OpWettingAngleRhs(const std::string field_name,
                    boost::shared_ptr<MatrixDouble> grad_h_ptr,
                    boost::shared_ptr<Range> ents_ptr = nullptr,
                    double wetting_angle = 0)
      : AssemblyBoundaryEleOp(field_name, field_name,
                              AssemblyBoundaryEleOp::OPROW),
        gradHPtr(grad_h_ptr), entsPtr(ents_ptr), wettingAngle(wetting_angle) {}
 
  MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data) {
    MoFEMFunctionBegin;
    if (entsPtr) {
      if (entsPtr->find(AssemblyBoundaryEleOp::getFEEntityHandle()) ==
          entsPtr->end())
        MoFEMFunctionReturnHot(0);
    }
    const double area = getMeasure();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_row_base = row_data.getFTensor0N();
    auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto s = wetting_angle_sub_stepping(getFEMethod()->ts_step);
 
    for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
      const double r = t_coords(0);
      const double alpha = t_w * cylindrical(r) * area;
      const double h_grad_norm = sqrt(t_grad_h(i) * t_grad_h(i) + eps); // ||∇h||
      const double cos_angle = std::cos(M_PI * wettingAngle / 180);
      const double rhs_wetting = s * eta2 * h_grad_norm * cos_angle; // s * η² * ||∇h|| * cos(angle)
 
      // cerr << "pass "
      //      << h_grad_norm <<"\n";
      int bb = 0;
      for (; bb != nbRows; ++bb) {
        locF[bb] += alpha * t_row_base * rhs_wetting;
        ++t_row_base;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb)
        ++t_row_base;
 
      ++t_w;
      ++t_grad_h;
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> gradHPtr;
  boost::shared_ptr<Range> entsPtr;
  double wettingAngle;
};
//! [OpWettingAngleRhs]
 
//! [OpNormalConstrainLhs]
/**
 * @brief Boundary normal constraint LHS (Jacobian block)
 * @param field_name_row Name of the field associated with the row operator
 * @param field_name_col Name of the field associated with the column operator
 * @return MoFEMErrorCode
 * 
 * Computes the matrix entries for ∫_Γ ​w_L​⋅(n⋅u) dΓ, and its transpose 
 * ∫_Γ ​w_u​⋅(λn) dΓ, where:
 * - w_L  : test function for Lagrange multiplier (λ)
 * - w_u  : test function for velocity (u)
 * - n    : outward normal vector on the boundary.
 */
struct OpNormalConstrainLhs : public AssemblyBoundaryEleOp {
 
  OpNormalConstrainLhs(const std::string field_name_row,
                       const std::string field_name_col)
      : AssemblyBoundaryEleOp(field_name_row, field_name_col,
                              AssemblyBoundaryEleOp::OPROWCOL) {
    assembleTranspose = true;
    sYmm = false;
  }
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_normal = getFTensor1Normal();
    auto t_row_base = row_data.getFTensor0N();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      auto t_mat = getFTensor1FromPtr<U_FIELD_DIM>(&locMat(0, 0));
 
      const double r = t_coords(0);
      const double alpha = t_w * cylindrical(r);
 
      int rr = 0;
      for (; rr != nbRows; ++rr) {
 
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        const double a = alpha * t_row_base;
 
        for (int cc = 0; cc != nbCols / U_FIELD_DIM; ++cc) {
          t_mat(i) += (a * t_col_base) * t_normal(i); // α * N_row * N_col * n
          ++t_col_base;
          ++t_mat;
        }
        ++t_row_base;
      }
 
      for (; rr < nbRowBaseFunctions; ++rr)
        ++t_row_base;
 
      ++t_w;
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  };
};
//! [OpNormalConstrainLhs]
 
//! [OpWettingAngleLhs]
/**
* @brief Boundary wetting-angle LHS (linearization)
* @param field_name Name of the field associated with the operator
* @param grad_h_ptr Pointer to the gradient of the free surface height matrix
* @param col_ind_ptr Pointer to the column indices vector
* @param col_diff_base_ptr Pointer to the column difference base functions vector
* @param ents_ptr Pointer to the range of entities (edges) to consider
* @param wetting_angle Contact wetting angle in degrees
* @return MoFEMErrorCode
* 
* Linearization of the wetting-angle boundary term (3.1d - Gamma part) used to assemble the
* corresponding Jacobian contributions for G (and coupling into H).
*
* δ R_g / δ h = ∫_Γ s ε² (∇h / ||∇h||) · ∇(δ h) cos(α) dΓ
*/
struct OpWettingAngleLhs : public BoundaryEleOp {
 
  OpWettingAngleLhs(
      const std::string row_field_name,
      boost::shared_ptr<MatrixDouble> grad_h_ptr,
      boost::shared_ptr<std::vector<VectorInt>> col_ind_ptr,
      boost::shared_ptr<std::vector<MatrixDouble>> col_diff_base_ptr,
      boost::shared_ptr<Range> ents_ptr = nullptr, double wetting_angle = 0)
      : BoundaryEleOp(row_field_name, BoundaryEleOp::OPROW),
        gradHPtr(grad_h_ptr), colIndicesPtr(col_ind_ptr),
        colDiffBaseFunctionsPtr(col_diff_base_ptr), entsPtr(ents_ptr),
        wettingAngle(wetting_angle) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) {
    MoFEMFunctionBegin;
    if (entsPtr) {
      if (entsPtr->find(BoundaryEleOp::getFEEntityHandle()) == entsPtr->end())
        MoFEMFunctionReturnHot(0);
    }
    const double area = getMeasure();
 
    const auto row_size = data.getIndices().size();
    if (row_size == 0)
      MoFEMFunctionReturnHot(0);
 
    auto integrate = [&](auto col_indicies, auto &col_diff_base_functions) {
      MoFEMFunctionBegin;
 
      const auto col_size = col_indicies.size();
 
      locMat.resize(row_size, col_size, false);
      locMat.clear();
      int nb_gp = getGaussPts().size2();
      int nb_rows = data.getIndices().size();
 
      auto t_w = getFTensor0IntegrationWeight();
      auto t_coords = getFTensor1CoordsAtGaussPts();
      auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
      auto t_row_base = data.getFTensor0N();
      int nb_row_base_functions = data.getN().size2();
 
      auto s = wetting_angle_sub_stepping(getFEMethod()->ts_step);
 
      for (int gg = 0; gg != nb_gp; ++gg) {
 
        const double r = t_coords(0);
        const double alpha = t_w * area * cylindrical(r);
        const double h_grad_norm = sqrt(t_grad_h(i) * t_grad_h(i) + eps);
        const double one_over_h_grad_norm = 1. / h_grad_norm;
        const double beta = s * alpha * eta2 * one_over_h_grad_norm *
                            std::cos(M_PI * wettingAngle / 180); // s * α * η² * (1/||∇h||) * cos(angle)
 
        int rr = 0;
        for (; rr != nb_rows; ++rr) {
          const double delta = beta * t_row_base; // β * N_row
 
          auto ptr = &col_diff_base_functions(gg, 0);
          auto t_col_diff_base = getFTensor1FromPtr<SPACE_DIM>(ptr);
 
          for (int cc = 0; cc != col_size; ++cc) {
            locMat(rr, cc) += t_col_diff_base(i) * (delta * t_grad_h(i));
            ++t_col_diff_base;
          }
          ++t_row_base;
        }
 
        for (; rr < nb_row_base_functions; ++rr) {
          ++t_row_base;
        }
 
        ++t_grad_h;
        ++t_w;
        ++t_coords;
      }
 
      MoFEMFunctionReturn(0);
    };
 
    for (auto c = 0; c != colIndicesPtr->size(); ++c) {
 
      auto &col_ind = (*colIndicesPtr)[c];
      if (col_ind.size()) {
        auto &diff_base = (*colDiffBaseFunctionsPtr)[c];
 
        CHKERR integrate(col_ind, diff_base);
 
        CHKERR MatSetValues(getKSPB(), data.getIndices().size(),
                            &*data.getIndices().begin(), col_ind.size(),
                            &*col_ind.begin(), &locMat(0, 0), ADD_VALUES);
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  MatrixDouble locMat;
 
  boost::shared_ptr<MatrixDouble> gradHPtr;
  boost::shared_ptr<Range> entsPtr;
  double wettingAngle;
  boost::shared_ptr<std::vector<VectorInt>> colIndicesPtr;
  boost::shared_ptr<std::vector<MatrixDouble>> colDiffBaseFunctionsPtr;
};
//! [OpWettingAngleLhs]
 
//! [OpRhsU]
/**
 * @brief Rhs for U (momentum residual)
 * @param field_name Name of the field associated with the operator
 * @param dot_u_ptr Pointer to the velocity time derivative matrix (∂U/∂t)
 * @param u_ptr Pointer to the velocity matrix (U)
 * @param grad_u_ptr Pointer to the velocity gradient matrix (∇U)
 * @param h_ptr Pointer to the free surface height vector (h)
 * @param grad_h_ptr Pointer to the free surface height gradient matrix (∇h)
 * @param g_ptr Pointer to the chemical potential vector (g)
 * @param p_ptr Pointer to the pressure vector (p)
 * @return MoFEMErrorCode
 *
 * This operator implements the volume terms of the momentum equation (not Jacobian)
 * (Lovric 3.1a)
 * 
 * R_u = ∫_Ω w·(ρ(h)(∂u/∂t + u·∇u + a_0) - κg∇h) - (∇·w)p + (∇w):(2μ(h)D(u))dΩ
 *
 * Code ↔ PDE term mapping (paper notation : code variable)
 * - ρ(ϕ) ∂u/∂t           : t_inertia_force(i) = (rho * alpha) * t_dot_u(i)
 * - ρ(ϕ) (u·∇)u         : t_convection(i)  = (rho * alpha) * (t_u(j) * t_grad_u(i,j))
 * - ∇·(2 μ(ϕ) D(u))      : t_stress built from get_D(2*mu) and t_grad_u,
 *                          then assembled with t_diff_base into locF
 * - −∇p (pressure grad)  : incorporated via t_kd(i,j)*t_p term inside t_stress
 * - −κ η ∇ϕ (surface)    : t_phase_force(i) = -alpha * kappa * t_g * t_grad_h(i)
 * - body force ρ a0      : t_gravity(SPACE_DIM-1) = alpha * rho * a0
 * - buoyancy force        : (currently commented out) t_buoyancy(SPACE_DIM-1) = -(alpha * rho * a0) * t_h;
 */
struct OpRhsU : public AssemblyDomainEleOp {
 
  OpRhsU(const std::string field_name,
         boost::shared_ptr<MatrixDouble> dot_u_ptr, // ∂U/∂t
         boost::shared_ptr<MatrixDouble> u_ptr,
         boost::shared_ptr<MatrixDouble> grad_u_ptr, // ∇U
         boost::shared_ptr<VectorDouble> h_ptr,
         boost::shared_ptr<MatrixDouble> grad_h_ptr, // ∇h
         boost::shared_ptr<VectorDouble> g_ptr, // chemical potential
         boost::shared_ptr<VectorDouble> p_ptr) // pressure
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        dotUPtr(dot_u_ptr), uPtr(u_ptr), gradUPtr(grad_u_ptr), hPtr(h_ptr),
        gradHPtr(grad_h_ptr), gPtr(g_ptr), pPtr(p_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_dot_u = getFTensor1FromMat<U_FIELD_DIM>(*dotUPtr);
    auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
    auto t_p = getFTensor0FromVec(*pPtr);
    auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*gradUPtr);
    auto t_h = getFTensor0FromVec(*hPtr);
    auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
    auto t_g = getFTensor0FromVec(*gPtr);
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto t_base = data.getFTensor0N();
    auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
 
    auto t_w = getFTensor0IntegrationWeight();
 
    constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<double>();
    FTensor::Tensor2_symmetric<double, SPACE_DIM> t_stress;
    FTensor::Tensor1<double, U_FIELD_DIM> t_phase_force;
    FTensor::Tensor1<double, U_FIELD_DIM> t_inertia_force;
    FTensor::Tensor1<double, U_FIELD_DIM> t_convection;
    FTensor::Tensor1<double, U_FIELD_DIM> t_buoyancy;
    FTensor::Tensor1<double, U_FIELD_DIM> t_gravity;
    FTensor::Tensor1<double, U_FIELD_DIM> t_forces;
 
    t_buoyancy(i) = 0;
    t_gravity(i) = 0;
 
    for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
      const double r = t_coords(0);
      const double alpha = t_w * vol * cylindrical(r);
 
      const double rho = phase_function(t_h, rho_diff, rho_ave);
      const double mu = phase_function(t_h, mu_diff, mu_ave);
 
      auto t_D = get_D(2 * mu);
 
      t_inertia_force(i) = (rho * alpha) * (t_dot_u(i)); // ρ(h)∂U/∂t
      // t_buoyancy(SPACE_DIM - 1) = -(alpha * rho * a0) * t_h;
      t_gravity(SPACE_DIM - 1) = (alpha * rho * a0); // ρ(h)a_0 (a_0 is gravity)
      t_phase_force(i) = -alpha * kappa * t_g * t_grad_h(i); // -κg∇h
      t_convection(i) = (rho * alpha) * (t_u(j) * t_grad_u(i, j)); // ρ(h)(U·∇)U
 
      t_stress(i, j) =
          alpha * (t_D(i, j, k, l) * t_grad_u(k, l) + t_kd(i, j) * t_p); // ∇·(2μ(h)D(U))
 
      auto t_nf = getFTensor1FromArray<U_FIELD_DIM, U_FIELD_DIM>(locF);
 
      t_forces(i) = t_inertia_force(i) + t_buoyancy(i) + t_gravity(i) +
                    t_convection(i) + t_phase_force(i);
 
      int bb = 0;
      for (; bb != nbRows / U_FIELD_DIM; ++bb) {
 
        // Assemble total force contributions into local residual
        t_nf(i) += t_base * t_forces(i);
        t_nf(i) += t_diff_base(j) * t_stress(i, j); // ∇w : ∇·(2μ(h)D(U)) - ∇p
 
        if (coord_type == CYLINDRICAL) {
          t_nf(0) += (t_base * (alpha / t_coords(0))) * (2 * mu * t_u(0) + t_p);
        }
 
        ++t_base;
        ++t_diff_base;
        ++t_nf;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb) {
        ++t_diff_base;
        ++t_base;
      }
 
      ++t_dot_u;
      ++t_u;
      ++t_grad_u;
      ++t_h;
      ++t_grad_h;
      ++t_g;
      ++t_p;
 
      ++t_w;
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> dotUPtr;
  boost::shared_ptr<MatrixDouble> uPtr;
  boost::shared_ptr<MatrixDouble> gradUPtr;
  boost::shared_ptr<VectorDouble> hPtr;
  boost::shared_ptr<MatrixDouble> gradHPtr;
  boost::shared_ptr<VectorDouble> gPtr;
  boost::shared_ptr<VectorDouble> pPtr;
};
//! [OpRhsU]
 
//! [OpLhsU_dU]
/**
 * @brief Lhs for U dU (Jacobian block ∂R_U/∂U)
 * @param field_name Name of the field associated with the operator
 * @param u_ptr Pointer to the velocity matrix (U)
 * @param grad_u_ptr Pointer to the velocity gradient matrix (∇U)
 * @param h_ptr Pointer to the free surface height vector (h)
 * @return MoFEMErrorCode
 *
 * Implements linearization of momentum residual w.r.t. velocity U:
 *
 * Matches the Jacobian terms required to solve the implicit momentum equation
 * used when solving the coupled NS–CH system (Lovric 3.1a linearization).
 */
struct OpLhsU_dU : public AssemblyDomainEleOp {
 
  OpLhsU_dU(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
            boost::shared_ptr<MatrixDouble> grad_u_ptr,
            boost::shared_ptr<VectorDouble> h_ptr)
      : AssemblyDomainEleOp(field_name, field_name,
                            AssemblyDomainEleOp::OPROWCOL),
        uPtr(u_ptr), gradUPtr(grad_u_ptr), hPtr(h_ptr) {
    sYmm = false;
    assembleTranspose = false;
  }
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
    auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*gradUPtr);
    auto t_h = getFTensor0FromVec(*hPtr);
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto t_row_base = row_data.getFTensor0N();
    auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
 
    auto t_w = getFTensor0IntegrationWeight();
 
    auto get_mat = [&](const int rr) {
      return getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(locMat, rr);
    };
 
    auto ts_a = getTSa();
    constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<double>(); // t_kd(i,j)=1 iff i=j
 
    for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
      const double r = t_coords(0);
      const double alpha = t_w * vol * cylindrical(r);
      const double rho = phase_function(t_h, rho_diff, rho_ave);
      const double mu = phase_function(t_h, mu_diff, mu_ave);
 
      const double beta0 = alpha * rho;
      const double beta1 = beta0 * ts_a;
      auto t_D = get_D(2 * mu); // Fourth-order viscosity tensor
 
      int rr = 0;
      for (; rr != nbRows / U_FIELD_DIM; ++rr) {
 
        auto t_mat = get_mat(rr * U_FIELD_DIM);
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
        
        // Symmetric on last two indices
        // ie Christof(i,j,k) == Christof(i,k,j)
        FTensor::Christof<double, SPACE_DIM, SPACE_DIM> t_d_stress;
        t_d_stress(l, j, k) = t_D(i, j, k, l) * (alpha * t_row_diff_base(i));
 
        for (int cc = 0; cc != nbCols / U_FIELD_DIM; ++cc) {// Assemble contributions for each col basis function
 
          const double bb = t_row_base * t_col_base;
 
          t_mat(i, j) += (beta1 * bb) * t_kd(i, j); // Derivative of inertia term ρ ∂U/∂t
          t_mat(i, j) += (beta0 * bb) * t_grad_u(i, j); // Derivative of convection term ρ (U·∇)U
          t_mat(i, j) +=
              (beta0 * t_row_base) * t_kd(i, j) * (t_col_diff_base(k) * t_u(k));
          t_mat(i, j) += t_d_stress(i, j, k) * t_col_diff_base(k); // Derivative of viscous term ∇·(2μD(U))
 
          if (coord_type == CYLINDRICAL) {
            t_mat(0, 0) += (bb * (alpha / t_coords(0))) * (2 * mu);
          }
 
          ++t_mat;
          ++t_col_base;
          ++t_col_diff_base;
        }
 
        ++t_row_base;
        ++t_row_diff_base;
      }
 
      for (; rr < nbRowBaseFunctions; ++rr) {
        ++t_row_diff_base;
        ++t_row_base;
      }
 
      ++t_u;
      ++t_grad_u;
      ++t_h;
 
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> uPtr;
  boost::shared_ptr<MatrixDouble> gradUPtr;
  boost::shared_ptr<VectorDouble> hPtr;
};
//! [OpLhsU_dU]
 
struct OpLoopSideGetDataForSideEle : ForcesAndSourcesCore::UserDataOperator {
 
  using UDO = ForcesAndSourcesCore::UserDataOperator;
 
  OpLoopSideGetDataForSideEle(
      const std::string field_name,
      boost::shared_ptr<std::vector<VectorInt>> col_indices_ptr,
      boost::shared_ptr<std::vector<MatrixDouble>> col_diff_basefunctions_ptr)
      : UDO(field_name, UDO::OPCOL), colIndicesPtr(col_indices_ptr),
        colDiffBaseFunctionsPtr(col_diff_basefunctions_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        DataForcesAndSourcesCore::EntData &data) {
    MoFEMFunctionBegin;
 
    if (type == MBVERTEX) {
      colIndicesPtr->clear();
      colDiffBaseFunctionsPtr->clear();
    }
 
    colIndicesPtr->push_back(data.getIndices());
    colDiffBaseFunctionsPtr->push_back(data.getDiffN());
 
    MoFEMFunctionReturn(0);
  }
 
protected:
  boost::shared_ptr<std::vector<VectorInt>> colIndicesPtr;
  boost::shared_ptr<std::vector<MatrixDouble>> colDiffBaseFunctionsPtr;
};
 
/**
 * @brief Lhs for U dH (Jacobian block ∂R_U/∂H)
 * @param field_name_u Name of the field associated with the row operator (U)
 * @param field_name_h Name of the field associated with the column operator (H)
 * @param dot_u_ptr Pointer to the velocity time derivative matrix (∂U/∂t)
 * @param u_ptr Pointer to the velocity matrix (U)
 * @param grad_u_ptr Pointer to the velocity gradient matrix (∇U)
 * @param h_ptr Pointer to the free surface height vector (h)
 * @param g_ptr Pointer to the chemical potential vector (g)
 * @return MoFEMErrorCode
 *
 * This operator builds the sensitivity of the momentum residual to changes
 * in the phase field h (paper 3.2a cross-coupling). It implements:
 *
 * In code: rho_dh = d_phase_function_h(t_h, rho_diff), mu_dh = d_phase_function_h(t_h, mu_diff)
 */
struct OpLhsU_dH : public AssemblyDomainEleOp {
 
  OpLhsU_dH(const std::string field_name_u, const std::string field_name_h,
            boost::shared_ptr<MatrixDouble> dot_u_ptr,
            boost::shared_ptr<MatrixDouble> u_ptr,
            boost::shared_ptr<MatrixDouble> grad_u_ptr,
            boost::shared_ptr<VectorDouble> h_ptr,
            boost::shared_ptr<VectorDouble> g_ptr)
      : AssemblyDomainEleOp(field_name_u, field_name_h,
                            AssemblyDomainEleOp::OPROWCOL),
        dotUPtr(dot_u_ptr), uPtr(u_ptr), gradUPtr(grad_u_ptr), hPtr(h_ptr),
        gPtr(g_ptr) {
    sYmm = false;
    assembleTranspose = false;
  }
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_dot_u = getFTensor1FromMat<U_FIELD_DIM>(*dotUPtr);
    auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
    auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*gradUPtr);
    auto t_h = getFTensor0FromVec(*hPtr);
    auto t_g = getFTensor0FromVec(*gPtr);
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto t_row_base = row_data.getFTensor0N();
    auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
 
    auto t_w = getFTensor0IntegrationWeight();
 
    FTensor::Tensor2_symmetric<double, SPACE_DIM> t_stress_dh;
    FTensor::Tensor1<double, U_FIELD_DIM> t_phase_force_dh;
    FTensor::Tensor1<double, U_FIELD_DIM> t_inertia_force_dh;
    FTensor::Tensor1<double, U_FIELD_DIM> t_convection_dh;
    FTensor::Tensor1<double, U_FIELD_DIM> t_buoyancy_dh;
    FTensor::Tensor1<double, U_FIELD_DIM> t_gravity_dh;
    FTensor::Tensor1<double, U_FIELD_DIM> t_forces_dh;
 
    t_buoyancy_dh(i) = 0;
    t_gravity_dh(i) = 0;
 
    for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
      const double r = t_coords(0);
      const double alpha = t_w * vol * cylindrical(r);
 
      const double rho_dh = d_phase_function_h(t_h, rho_diff); // dρ/dh
      const double mu_dh = d_phase_function_h(t_h, mu_diff); // dμ/dh
 
      auto t_D_dh = get_D(2 * mu_dh);
 
      t_inertia_force_dh(i) = (alpha * rho_dh) * t_dot_u(i); // dρ/dh * ∂U/∂t
      // t_buoyancy_dh(SPACE_DIM - 1) = -(alpha * a0) * (rho + rho_dh * t_h);
      t_gravity_dh(SPACE_DIM - 1) = (alpha * rho_dh * a0); // dρ/dh * a_0
      t_convection_dh(i) = (rho_dh * alpha) * (t_u(j) * t_grad_u(i, j)); // dρ/dh * (U·∇)U
      const double t_phase_force_g_dh = -alpha * kappa * t_g; // -κ g
      t_forces_dh(i) = t_inertia_force_dh(i) + t_buoyancy_dh(i) +
                       t_gravity_dh(i) + t_convection_dh(i); // Total dR_U/dh force
 
      t_stress_dh(i, j) = alpha * (t_D_dh(i, j, k, l) * t_grad_u(k, l)); // d(∇·(2μD(U)))/dh
 
      int rr = 0;
      for (; rr != nbRows / U_FIELD_DIM; ++rr) {
 
        auto t_mat =
            getFTensor1FromMat<U_FIELD_DIM, 1>(locMat, rr * U_FIELD_DIM);
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
 
        for (int cc = 0; cc != nbCols; ++cc) {
 
          const double bb = t_row_base * t_col_base;
          t_mat(i) += t_forces_dh(i) * bb; // Assemble total dR_U/dh force contribution
          t_mat(i) += (t_phase_force_g_dh * t_row_base) * t_col_diff_base(i); // -κ g ∇h term
          t_mat(i) += (t_row_diff_base(j) * t_col_base) * t_stress_dh(i, j); // Viscous term contribution via dμ/dh
 
          if (coord_type == CYLINDRICAL) {
            t_mat(0) += (bb * (alpha / t_coords(0))) * (2 * mu_dh * t_u(0));
          }
 
          ++t_mat;
          ++t_col_base;
          ++t_col_diff_base;
        }
 
        ++t_row_base;
        ++t_row_diff_base;
      }
 
      for (; rr < nbRowBaseFunctions; ++rr) {
        ++t_row_diff_base;
        ++t_row_base;
      }
 
      ++t_dot_u;
      ++t_u;
      ++t_grad_u;
      ++t_h;
      ++t_g;
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> dotUPtr;
  boost::shared_ptr<MatrixDouble> uPtr;
  boost::shared_ptr<MatrixDouble> gradUPtr;
  boost::shared_ptr<VectorDouble> hPtr;
  boost::shared_ptr<VectorDouble> gPtr;
};
 
/**
 * @brief Lhs for U dG (Jacobian block ∂R_U/∂G)
 * @param field_name_u Name of the field associated with the row operator (U)
 * @param field_name_h Name of the field associated with the column operator (H)
 * @param grad_h_ptr Pointer to the free surface height gradient matrix (∇h)
 * @return MoFEMErrorCode
 *
 * This operator assembles the contribution of the chemical potential G into
 * the momentum Jacobian. It corresponds to linearization of the surface
 * tension term −κ g ∇h wrt g.
 */
struct OpLhsU_dG : public AssemblyDomainEleOp {
 
  OpLhsU_dG(const std::string field_name_u, const std::string field_name_h,
            boost::shared_ptr<MatrixDouble> grad_h_ptr)
      : AssemblyDomainEleOp(field_name_u, field_name_h,
                            AssemblyDomainEleOp::OPROWCOL),
        gradHPtr(grad_h_ptr) {
    sYmm = false;
    assembleTranspose = false;
  }
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto t_row_base = row_data.getFTensor0N();
    auto t_w = getFTensor0IntegrationWeight();
 
    for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
      const double r = t_coords(0);
      const double alpha = t_w * vol * cylindrical(r);
 
      FTensor::Tensor1<double, SPACE_DIM> t_phase_force_dg;
      t_phase_force_dg(i) = -alpha * kappa * t_grad_h(i); // -κ ∇h
 
      int rr = 0;
      for (; rr != nbRows / U_FIELD_DIM; ++rr) {
        auto t_mat =
            getFTensor1FromMat<U_FIELD_DIM, 1>(locMat, rr * U_FIELD_DIM);
        auto t_col_base = col_data.getFTensor0N(gg, 0);
 
        for (int cc = 0; cc != nbCols; ++cc) {
          const double bb = t_row_base * t_col_base;
          t_mat(i) += t_phase_force_dg(i) * bb;
 
          ++t_mat;
          ++t_col_base;
        }
 
        ++t_row_base;
      }
 
      for (; rr < nbRowBaseFunctions; ++rr) {
        ++t_row_base;
      }
 
      ++t_grad_h;
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> gradHPtr;
};
 
/**
 * @brief Rhs for H (phase-field residual)
 * @param field_name Name of the field associated with the operator
 * @param u_ptr Pointer to the velocity matrix (U)
 * @param dot_h_ptr Pointer to the free surface height time derivative vector (∂H/∂t)
 * @param h_ptr Pointer to the free surface height vector (H)
 * @param grad_h_ptr Pointer to the free surface height gradient matrix (∇H)
 * @param grad_g_ptr Pointer to the chemical potential gradient matrix (∇G)
 * @return MoFEMErrorCode
 *
 * Maps the conserved phase evolution (Lovric 3.1c):
 * R_h = ∫_Ω v(∂h/∂t + u·∇h) + ∇v·(M(h) ∇g)) dΩ
 *
 * The template boolean I selects initialization (I=true) vs evolution (I=false).
 */
template <bool I> struct OpRhsH : public AssemblyDomainEleOp {
 
  OpRhsH(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
         boost::shared_ptr<VectorDouble> dot_h_ptr,
         boost::shared_ptr<VectorDouble> h_ptr,
         boost::shared_ptr<MatrixDouble> grad_h_ptr,
         boost::shared_ptr<MatrixDouble> grad_g_ptr)
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        uPtr(u_ptr), dotHPtr(dot_h_ptr), hPtr(h_ptr), gradHPtr(grad_h_ptr),
        gradGPtr(grad_g_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_base = data.getFTensor0N();
    auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
 
#ifndef NDEBUG
    if (data.getDiffN().size1() != data.getN().size1())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
    if (data.getDiffN().size2() != data.getN().size2() * SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::error)
          << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
      MOFEM_LOG("SELF", Sev::error) << data.getN();
      MOFEM_LOG("SELF", Sev::error) << data.getDiffN();
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
    }
#endif
 
    if constexpr (I) {
      // Find h ∈ H¹(Ω) to minimise ∫(h - h_init)² dΩ + ∫∇v·(M ∇g) dΩ
 
      auto t_h = getFTensor0FromVec(*hPtr);
      auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
 
      for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
        const double r = t_coords(0);
        const double alpha = t_w * vol * cylindrical(r);
 
        const double set_h = init_h(t_coords(0), t_coords(1), t_coords(2));
        const double m = get_M(set_h) * alpha;
 
        int bb = 0;
        for (; bb != nbRows; ++bb) {
          locF[bb] += (t_base * alpha) * (t_h - set_h); // minimize ∫(h - h_init)² dΩ
          locF[bb] += (t_diff_base(i) * m) * t_grad_g(i);
          ++t_base;
          ++t_diff_base;
        }
 
        for (; bb < nbRowBaseFunctions; ++bb) {
          ++t_base;
          ++t_diff_base;
        }
 
        ++t_h;
        ++t_grad_g;
 
        ++t_coords;
        ++t_w;
      }
 
    } else {
 
      auto t_dot_h = getFTensor0FromVec(*dotHPtr);
      auto t_h = getFTensor0FromVec(*hPtr);
      auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
      auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
      auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
 
      for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
        const double r = t_coords(0);
        const double alpha = t_w * vol * cylindrical(r);
 
        const double m = get_M(t_h) * alpha;
 
        int bb = 0;
        for (; bb != nbRows; ++bb) {
          locF[bb] += (t_base * alpha) * (t_dot_h); // time derivative ∂h/∂t
          locF[bb] += (t_base * alpha) * (t_grad_h(i) * t_u(i)); // convection u·∇h
          locF[bb] += (t_diff_base(i) * t_grad_g(i)) * m; // mobility/divergence (M ∇g)
          ++t_base;
          ++t_diff_base;
        }
 
        for (; bb < nbRowBaseFunctions; ++bb) {
          ++t_base;
          ++t_diff_base;
        }
 
        ++t_dot_h;
        ++t_h;
        ++t_grad_g;
        ++t_u;
        ++t_grad_h;
 
        ++t_coords;
        ++t_w;
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> uPtr;
  boost::shared_ptr<VectorDouble> dotHPtr;
  boost::shared_ptr<VectorDouble> hPtr;
  boost::shared_ptr<MatrixDouble> gradHPtr;
  boost::shared_ptr<MatrixDouble> gradGPtr;
};
 
/**
 * @brief Lhs for H dU
 * @param h_field_name Name of the field associated with the row operator (H)
 * @param u_field_name Name of the field associated with the column operator (U)
 * @param grad_h_ptr Pointer to the free surface height gradient matrix (∇H)
 * @return MoFEMErrorCode
 * 
 * (Jacobian block ∂R_H/∂U)
 * Linearization of the conserved phase equation (Lovric 3.1c) w.r.t. U
 */
struct OpLhsH_dU : public AssemblyDomainEleOp {
  OpLhsH_dU(const std::string h_field_name, const std::string u_field_name,
            boost::shared_ptr<MatrixDouble> grad_h_ptr)
      : AssemblyDomainEleOp(h_field_name, u_field_name,
                            AssemblyDomainEleOp::OPROWCOL),
        gradHPtr(grad_h_ptr) {
    sYmm = false;
    assembleTranspose = false;
  }
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto t_row_base = row_data.getFTensor0N();
    auto t_w = getFTensor0IntegrationWeight();
 
    for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
      const auto r = t_coords(0);
      const auto alpha = t_w * vol * cylindrical(r);
      auto t_mat = getFTensor1FromPtr<U_FIELD_DIM>(&locMat(0, 0));
      FTensor::Tensor1<double, SPACE_DIM> t_row;
 
      int rr = 0;
      for (; rr != nbRows; ++rr) {
        t_row(i) = (alpha * t_row_base) * t_grad_h(i);
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        for (int cc = 0; cc != nbCols / U_FIELD_DIM; ++cc) {
          t_mat(i) += t_row(i) * t_col_base;
          ++t_mat;
          ++t_col_base;
        }
        ++t_row_base;
      }
 
      for (; rr < nbRowBaseFunctions; ++rr)
        ++t_row_base;
 
      ++t_grad_h;
      ++t_w;
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> gradHPtr;
};
 
/**
 * @brief Lhs for H dH (Jacobian block ∂R_H/∂H)
 * @param field_name Name of the field associated with the operator
 * @param u_ptr Pointer to the velocity matrix (U)
 * @param h_ptr Pointer to the free surface height vector (H)
 * @param grad_g_ptr Pointer to the chemical potential gradient matrix (∇G)
 * @return MoFEMErrorCode
 *
 * Linearization of the conserved phase equation (3.1c) w.r.t. H
 * - time-stepping mass scaling (ts_a) terms
 * - contributions from mobility derivative M'(h)
 * - coupling terms with chemical potential gradient where M depends on h
 *
 * When I==true the operator builds initialization-mode blocks (uses init_h)
 */
template <bool I> struct OpLhsH_dH : public AssemblyDomainEleOp {
 
  OpLhsH_dH(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
            boost::shared_ptr<VectorDouble> h_ptr,
            boost::shared_ptr<MatrixDouble> grad_g_ptr)
      : AssemblyDomainEleOp(field_name, field_name,
                            AssemblyDomainEleOp::OPROWCOL),
        uPtr(u_ptr), hPtr(h_ptr), gradGPtr(grad_g_ptr) {
    sYmm = false;
  }
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
    auto t_row_base = row_data.getFTensor0N();
    auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
 
#ifndef NDEBUG
    if (row_data.getDiffN().size1() != row_data.getN().size1())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
    if (row_data.getDiffN().size2() != row_data.getN().size2() * SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::error)
          << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
      MOFEM_LOG("SELF", Sev::error) << row_data.getN();
      MOFEM_LOG("SELF", Sev::error) << row_data.getDiffN();
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
    }
 
    if (col_data.getDiffN().size1() != col_data.getN().size1())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
    if (col_data.getDiffN().size2() != col_data.getN().size2() * SPACE_DIM) {
      MOFEM_LOG("SELF", Sev::error)
          << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
      MOFEM_LOG("SELF", Sev::error) << col_data.getN();
      MOFEM_LOG("SELF", Sev::error) << col_data.getDiffN();
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
    }
#endif
 
    if constexpr (I) {
 
      auto t_h = getFTensor0FromVec(*hPtr);
      auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
 
      for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
        const double r = t_coords(0);
        const double alpha = t_w * vol * cylindrical(r);
 
        int rr = 0;
        for (; rr != nbRows; ++rr) {
 
          auto t_col_base = col_data.getFTensor0N(gg, 0);
          auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
 
          for (int cc = 0; cc != nbCols; ++cc) {
 
            locMat(rr, cc) += (t_row_base * t_col_base * alpha);
 
            ++t_col_base;
            ++t_col_diff_base;
          }
 
          ++t_row_base;
          ++t_row_diff_base;
        }
 
        for (; rr < nbRowBaseFunctions; ++rr) {
          ++t_row_base;
          ++t_row_diff_base;
        }
 
        ++t_h;
        ++t_grad_g;
        ++t_w;
        ++t_coords;
      }
 
    } else {
 
      auto t_h = getFTensor0FromVec(*hPtr);
      auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
      auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
 
      auto ts_a = getTSa();
 
      for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
        const double r = t_coords(0);
        const double alpha = t_w * vol * cylindrical(r);
 
        auto m_dh = get_M_dh(t_h) * alpha;
        int rr = 0;
        for (; rr != nbRows; ++rr) {
 
          auto t_col_base = col_data.getFTensor0N(gg, 0);
          auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
 
          for (int cc = 0; cc != nbCols; ++cc) {
 
            locMat(rr, cc) += (t_row_base * t_col_base * alpha) * ts_a; // time-stepping mass term
            locMat(rr, cc) +=
                (t_row_base * alpha) * (t_col_diff_base(i) * t_u(i)); // convection term u·∇h 
            locMat(rr, cc) +=
                (t_row_diff_base(i) * t_grad_g(i)) * (t_col_base * m_dh); // mobility derivative term M'(h) ∇g
 
            ++t_col_base;
            ++t_col_diff_base;
          }
 
          ++t_row_base;
          ++t_row_diff_base;
        }
 
        for (; rr < nbRowBaseFunctions; ++rr) {
          ++t_row_base;
          ++t_row_diff_base;
        }
 
        ++t_u;
        ++t_h;
        ++t_grad_g;
        ++t_w;
        ++t_coords;
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> uPtr;
  boost::shared_ptr<VectorDouble> hPtr;
  boost::shared_ptr<MatrixDouble> gradGPtr;
};
 
/**
 * @brief Lhs for H dG
 * @param field_name_h Name of the field associated with the row operator (H)
 * @param field_name_g Name of the field associated with the column operator (G)
 * @param h_ptr Pointer to the free surface height vector (H)
 * @return MoFEMErrorCode
 * 
 * (Jacobian block ∂R_H/∂G)
 * Linearization of the conserved phase equation (Lovric 3.1c) w.r.t. G
 */
template <bool I> struct OpLhsH_dG : public AssemblyDomainEleOp {
 
  OpLhsH_dG(const std::string field_name_h, const std::string field_name_g,
            boost::shared_ptr<VectorDouble> h_ptr)
      : AssemblyDomainEleOp(field_name_h, field_name_g,
                            AssemblyDomainEleOp::OPROWCOL),
        hPtr(h_ptr) {
    sYmm = false;
    assembleTranspose = false;
  }
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_h = getFTensor0FromVec(*hPtr);
 
    auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
      const double r = t_coords(0);
      const double alpha = t_w * vol * cylindrical(r);
 
      double set_h;
      if constexpr (I)
        set_h = init_h(t_coords(0), t_coords(1), t_coords(2));
      else
        set_h = t_h;
 
      auto m = get_M(set_h) * alpha;
 
      int rr = 0;
      for (; rr != nbRows; ++rr) {
        auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
 
        for (int cc = 0; cc != nbCols; ++cc) {
          locMat(rr, cc) += (t_row_diff_base(i) * t_col_diff_base(i)) * m; // mobility/divergence (M ∇g)
 
          ++t_col_diff_base;
        }
 
        ++t_row_diff_base;
      }
 
      for (; rr < nbRowBaseFunctions; ++rr) {
        ++t_row_diff_base;
      }
 
      ++t_h;
      ++t_w;
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<VectorDouble> hPtr;
};
 
/**
 * @brief Rhs for G (chemical potential residual)
 * @param field_name Name of the field associated with the operator
 * @param h_ptr Pointer to the free surface height vector (H)
 * @param grad_h_ptr Pointer to the free surface height gradient matrix (∇H)
 * @param g_ptr Pointer to the chemical potential vector (G)
 * @return MoFEMErrorCode
 *
 * Implements Lovric 3.1d (excpet for wetting boundary term - see OpWettingAngleRhs/Lhs):
 * R_g = ∫_Ω s(g - f(h)) - ∇s·(ε²∇h) dΩ
 * 
 * The template boolean I selects initialization (I=true) vs evolution (I=false).
 */
template <bool I> struct OpRhsG : public AssemblyDomainEleOp {
 
  OpRhsG(const std::string field_name, boost::shared_ptr<VectorDouble> h_ptr,
         boost::shared_ptr<MatrixDouble> grad_h_ptr,
         boost::shared_ptr<VectorDouble> g_ptr)
      : AssemblyDomainEleOp(field_name, field_name, AssemblyDomainEleOp::OPROW),
        hPtr(h_ptr), gradHPtr(grad_h_ptr), gPtr(g_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_h = getFTensor0FromVec(*hPtr);
    auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
    auto t_g = getFTensor0FromVec(*gPtr);
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto t_base = data.getFTensor0N();
    auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
    auto t_w = getFTensor0IntegrationWeight();
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      const double r = t_coords(0);
      const double alpha = t_w * vol * cylindrical(r);
 
      double set_h;
      if constexpr (I)
        set_h = init_h(t_coords(0), t_coords(1), t_coords(2));
      else
        set_h = t_h;
 
      const double f = get_f(set_h);
 
      int bb = 0;
      for (; bb != nbRows; ++bb) {
        locF[bb] += (t_base * alpha) * (t_g - f); // Bulk term: (g - f(h))
        locF[bb] -= (t_diff_base(i) * (eta2 * alpha)) * t_grad_h(i); // Diffusion term: -ε²∇h
        ++t_base;
        ++t_diff_base;
      }
 
      for (; bb < nbRowBaseFunctions; ++bb) {
        ++t_base;
        ++t_diff_base;
      }
 
      ++t_h;
      ++t_grad_h;
      ++t_g;
 
      ++t_coords;
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<VectorDouble> hPtr;
  boost::shared_ptr<MatrixDouble> gradHPtr;
  boost::shared_ptr<VectorDouble> gPtr;
};
 
/**
 * @brief Lhs for G dH
 * @param field_name_g Name of the field associated with the row operator (G)
 * @param field_name_h Name of the field associated with the column operator (H)
 * @param h_ptr Pointer to the free surface height vector (H)
 * @return MoFEMErrorCode
 * 
 * (Jacobian block ∂R_G/∂H)
 * Linearization of the chemical potential equation (Lovric 3.1d) w.r.t. H
 */
template <bool I> struct OpLhsG_dH : public AssemblyDomainEleOp {
 
  OpLhsG_dH(const std::string field_name_g, const std::string field_name_h,
            boost::shared_ptr<VectorDouble> h_ptr)
      : AssemblyDomainEleOp(field_name_g, field_name_h,
                            AssemblyDomainEleOp::OPROWCOL),
        hPtr(h_ptr) {
    sYmm = false;
    assembleTranspose = false;
  }
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
    auto t_h = getFTensor0FromVec(*hPtr);
 
    auto t_row_base = row_data.getFTensor0N();
    auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      const double r = t_coords(0);
      const double alpha = t_w * vol * cylindrical(r);
 
      const double f_dh = get_f_dh(t_h) * alpha; // derivative of f(h)
      const double beta = eta2 * alpha; // coefficient for diffusion term -ε²
 
      int rr = 0;
      for (; rr != nbRows; ++rr) {
 
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
 
        for (int cc = 0; cc != nbCols; ++cc) {
 
          if constexpr (I == false)
            locMat(rr, cc) -= (t_row_base * t_col_base) * f_dh; // bulk term derivative (g - f(h))
          locMat(rr, cc) -= (t_row_diff_base(i) * beta) * t_col_diff_base(i); // diffusion term derivative -ε²∇h
 
          ++t_col_base;
          ++t_col_diff_base;
        }
 
        ++t_row_base;
        ++t_row_diff_base;
      }
 
      for (; rr < nbRowBaseFunctions; ++rr) {
        ++t_row_base;
        ++t_row_diff_base;
      }
 
      ++t_h;
      ++t_w;
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<VectorDouble> hPtr;
};
 
/**
 * @brief Lhs for G dG
 * @param field_name Name of the field associated with the operator
 * @return MoFEMErrorCode
 * 
 * (Jacobian block ∂R_G/∂G)
 * Linearization of the chemical potential equation (Lovric 3.1d) w.r.t. G
 */
struct OpLhsG_dG : public AssemblyDomainEleOp {
 
  OpLhsG_dG(const std::string field_name)
      : AssemblyDomainEleOp(field_name, field_name,
                            AssemblyDomainEleOp::OPROWCOL) {
    sYmm = true;
  }
 
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
 
    const double vol = getMeasure();
 
    auto t_row_base = row_data.getFTensor0N();
    auto t_w = getFTensor0IntegrationWeight();
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      const double r = t_coords(0);
      const double alpha = t_w * vol * cylindrical(r);
 
      int rr = 0;
      for (; rr != nbRows; ++rr) {
        auto t_col_base = col_data.getFTensor0N(gg, 0);
        const double beta = alpha * t_row_base;
        for (int cc = 0; cc != nbCols; ++cc) {
          locMat(rr, cc) += (t_col_base * beta);
          ++t_col_base;
        }
 
        ++t_row_base;
      }
 
      for (; rr < nbRowBaseFunctions; ++rr) {
        ++t_row_base;
      }
 
      ++t_w;
      ++t_coords;
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
};
 
} // namespace FreeSurfaceOps