affine_tet_stress_gram.cpp

Go to the documentation of this file.

/**
 * \file affine_tet_stress_gram.cpp
 *
 * \brief One-tetrahedron diagnostic for the EP stress block.
 */
 
constexpr int adolc_tag = 1;
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <limits>
#include <numeric>
#include <sstream>
#include <string>
#include <vector>
 
#include <MatOps.hpp>
#include <MatElastic.hpp>
#include <EshelbianPlasticity.hpp>
using namespace EshelbianPlasticity;
 
static char help[] =
    "Assemble EP P+bubble stress diagnostics on one affine tet.\n"
    "  -diag_matrix_type gram|a00\n"
    "  -diag_include_bubble true|false\n"
    "  -diag_rank_tol <real>\n"
    "  -diag_check_full_rank true|false\n"
    "  -diag_view_matrix true|false\n"
    "  -broken_hdiv_base demkowicz|ainsworth\n\n";
 
/// Selects which diagnostic matrix is assembled.
enum class DiagnosticMatrixType { GRAM, A00 };
 
/// Reads -diag_matrix_type from PETSc options.
static MoFEMErrorCode getDiagnosticMatrixType(DiagnosticMatrixType &type) {
  MoFEMFunctionBegin;
 
  char matrix_type[32] = "gram";
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-diag_matrix_type",
                               matrix_type, sizeof(matrix_type),
                               PETSC_NULLPTR);
  const std::string type_name(matrix_type);
 
  if (type_name == "gram") {
    type = DiagnosticMatrixType::GRAM;
  } else if (type_name == "a00") {
    type = DiagnosticMatrixType::A00;
  } else {
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "Unknown -diag_matrix_type '%s'; use 'gram' or 'a00'",
            matrix_type);
  }
 
  MoFEMFunctionReturn(0);
}
 
/// Returns the PETSc option name for a diagnostic matrix type.
static const char *getDiagnosticMatrixTypeName(DiagnosticMatrixType type) {
  switch (type) {
  case DiagnosticMatrixType::GRAM:
    return "gram";
  case DiagnosticMatrixType::A00:
    return "a00";
  }
  return "unknown";
}
 
/// Stores a short rank summary for the assembled matrix.
struct RankDiagnostic {
  size_t rowsNb = 0;
  size_t colsNb = 0;
  size_t maxRank = 0;
  size_t rank = 0;
  double frobeniusNorm = 0;
  double symmetryDefect = 0;
  double absRankTol = 0;
  double minNonzeroSingular = 0;
  double maxSingular = 0;
};
 
/// Creates one unit affine tetrahedron mesh.
static MoFEMErrorCode createAffineTet(moab::Interface &moab,
                                      EntityHandle &meshset) {
  MoFEMFunctionBegin;
 
  double tet_coords[] = {0,   0,   0,   1.0, 0,   0,
                         0,   1.0, 0,   0,   0,   1.0};
  EntityHandle nodes[4];
  for (int nn = 0; nn < 4; nn++) {
    CHKERR moab.create_vertex(&tet_coords[3 * nn], nodes[nn]);
  }
 
  EntityHandle tet;
  CHKERR moab.create_element(MBTET, nodes, 4, tet);
  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
 
  Range tets;
  tets.insert(tet);
  CHKERR moab.add_entities(meshset, tets);
  for (int dim = 0; dim != 3; ++dim) {
    Range adj;
    CHKERR moab.get_adjacencies(tets, dim, true, adj,
                                moab::Interface::UNION);
    CHKERR moab.add_entities(meshset, adj);
  }
 
  MoFEMFunctionReturn(0);
}
 
/// Copies a PETSc matrix into a dense uBLAS matrix.
static MoFEMErrorCode getDenseMatrix(Mat A, MatrixDouble &dense) {
  MoFEMFunctionBegin;
 
  PetscInt rows_nb, cols_nb;
  CHKERR MatGetSize(A, &rows_nb, &cols_nb);
  dense.resize(rows_nb, cols_nb, false);
  dense.clear();
 
  std::vector<PetscInt> rows(rows_nb), cols(cols_nb);
  std::iota(rows.begin(), rows.end(), 0);
  std::iota(cols.begin(), cols.end(), 0);
 
  std::vector<PetscScalar> values(rows_nb * cols_nb, 0);
  CHKERR MatGetValues(A, rows_nb, rows.data(), cols_nb, cols.data(),
                      values.data());
  for (PetscInt rr = 0; rr != rows_nb; ++rr) {
    for (PetscInt cc = 0; cc != cols_nb; ++cc) {
      dense(rr, cc) = PetscRealPart(values[rr * cols_nb + cc]);
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
/// Prints a dense matrix row by row.
static MoFEMErrorCode printDenseMatrix(const MatrixDouble &dense) {
  MoFEMFunctionBegin;
  for (size_t rr = 0; rr != dense.size1(); ++rr) {
    std::ostringstream line;
    line << std::scientific << std::setprecision(16);
    for (size_t cc = 0; cc != dense.size2(); ++cc) {
      line << (cc ? " " : "") << dense(rr, cc);
    }
    MOFEM_LOG("EP", Sev::inform) << line.str();
  }
  MoFEMFunctionReturn(0);
}
 
/// Computes rank from the smaller of A*A^T and A^T*A.
static MoFEMErrorCode computeRankDiagnostic(const MatrixDouble &dense,
                                            const double rel_rank_tol,
                                            RankDiagnostic &rank_data) {
  MoFEMFunctionBegin;
 
  rank_data.rowsNb = dense.size1();
  rank_data.colsNb = dense.size2();
  rank_data.maxRank = std::min(rank_data.rowsNb, rank_data.colsNb);
 
  for (size_t rr = 0; rr != rank_data.rowsNb; ++rr) {
    for (size_t cc = 0; cc != rank_data.colsNb; ++cc) {
      const double a = dense(rr, cc);
      rank_data.frobeniusNorm += a * a;
      if (rank_data.rowsNb == rank_data.colsNb) {
        const double d = a - dense(cc, rr);
        rank_data.symmetryDefect += d * d;
      }
    }
  }
  rank_data.frobeniusNorm = std::sqrt(rank_data.frobeniusNorm);
  rank_data.symmetryDefect = std::sqrt(rank_data.symmetryDefect);
 
  if (!rank_data.maxRank)
    MoFEMFunctionReturnHot(0);
 
  const bool use_row_gram = rank_data.rowsNb <= rank_data.colsNb;
  MatrixDouble gram(rank_data.maxRank, rank_data.maxRank);
  gram.clear();
 
  if (use_row_gram) {
    for (size_t ii = 0; ii != rank_data.rowsNb; ++ii) {
      for (size_t jj = ii; jj != rank_data.rowsNb; ++jj) {
        double value = 0;
        for (size_t cc = 0; cc != rank_data.colsNb; ++cc)
          value += dense(ii, cc) * dense(jj, cc);
        gram(ii, jj) = value;
        gram(jj, ii) = value;
      }
    }
  } else {
    for (size_t ii = 0; ii != rank_data.colsNb; ++ii) {
      for (size_t jj = ii; jj != rank_data.colsNb; ++jj) {
        double value = 0;
        for (size_t rr = 0; rr != rank_data.rowsNb; ++rr)
          value += dense(rr, ii) * dense(rr, jj);
        gram(ii, jj) = value;
        gram(jj, ii) = value;
      }
    }
  }
 
  VectorDouble eig(rank_data.maxRank);
  MatrixDouble eig_vec(rank_data.maxRank, rank_data.maxRank);
  CHKERR computeEigenValuesSymmetric(gram, eig, eig_vec);
 
  double max_eig = 0;
  for (size_t ii = 0; ii != rank_data.maxRank; ++ii)
    max_eig = std::max(max_eig, std::abs(eig(ii)));
  rank_data.maxSingular = std::sqrt(std::max(max_eig, 0.0));
  rank_data.absRankTol =
      rel_rank_tol * std::max(rank_data.maxSingular, 1.0);
 
  const double eig_abs_tol = rank_data.absRankTol * rank_data.absRankTol;
  double min_nonzero = std::numeric_limits<double>::max();
  for (size_t ii = 0; ii != rank_data.maxRank; ++ii) {
    const double a = std::abs(eig(ii));
    if (a > eig_abs_tol) {
      ++rank_data.rank;
      min_nonzero = std::min(min_nonzero, std::sqrt(a));
    }
  }
  rank_data.minNonzeroSingular = rank_data.rank ? min_nonzero : 0.0;
 
  MoFEMFunctionReturn(0);
}
 
/// Prints a compact rank summary for the assembled matrix.
static MoFEMErrorCode analyseDenseStressBlock(
    const MatrixDouble &dense, const DiagnosticMatrixType matrix_type,
    const PetscBool include_bubble, const double rel_rank_tol,
    const PetscBool check_full_rank, const PetscBool view_matrix) {
  MoFEMFunctionBegin;
 
  RankDiagnostic rank_data;
  CHKERR computeRankDiagnostic(dense, rel_rank_tol, rank_data);
 
  MOFEM_LOG("EP", Sev::inform)
      << "diagnostic matrix type: " << getDiagnosticMatrixTypeName(matrix_type);
  if (matrix_type == DiagnosticMatrixType::GRAM) {
    MOFEM_LOG("EP", Sev::inform)
        << "stress basis: " << (include_bubble ? "P+bubble" : "P");
  } else {
    MOFEM_LOG("EP", Sev::inform)
        << "A00 fields: u, bubble, P, omega, wL2";
  }
  MOFEM_LOG("EP", Sev::inform)
      << "matrix size: " << rank_data.rowsNb << " x " << rank_data.colsNb;
  MOFEM_LOG("EP", Sev::inform)
      << std::scientific << std::setprecision(16)
      << "frobenius norm: " << rank_data.frobeniusNorm;
  if (rank_data.rowsNb == rank_data.colsNb) {
    MOFEM_LOG("EP", Sev::inform)
        << std::scientific << std::setprecision(16)
        << "symmetry defect ||A-A^T||_F: " << rank_data.symmetryDefect;
  }
  MOFEM_LOG("EP", Sev::inform)
      << std::scientific << std::setprecision(16)
      << "relative rank tolerance: " << rel_rank_tol;
  MOFEM_LOG("EP", Sev::inform)
      << std::scientific << std::setprecision(16)
      << "absolute singular tolerance: " << rank_data.absRankTol;
  MOFEM_LOG("EP", Sev::inform)
      << std::scientific << std::setprecision(16)
      << "singular value range approx: [" << rank_data.minNonzeroSingular
      << ", " << rank_data.maxSingular << "]";
  MOFEM_LOG("EP", Sev::inform)
      << "numerical rank: " << rank_data.rank << " / " << rank_data.maxRank;
 
  if (view_matrix) {
    MOFEM_LOG("EP", Sev::inform) << "dense diagnostic matrix:";
    CHKERR printDenseMatrix(dense);
  }
 
  if (check_full_rank && rank_data.rank != rank_data.maxRank) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "%s diagnostic block is rank deficient: rank %zu / %zu",
            getDiagnosticMatrixTypeName(matrix_type), rank_data.rank,
            rank_data.maxRank);
  }
 
  MoFEMFunctionReturn(0);
}
 
/// Builds the affine tet problem and assembles the diagnostic matrix.
static MoFEMErrorCode assembleStressDiagnostic(
    const DiagnosticMatrixType matrix_type, const PetscBool include_bubble,
    const double rel_rank_tol, const PetscBool check_full_rank,
    const PetscBool view_matrix) {
  MoFEMFunctionBegin;
 
  int comm_size = 0;
  MPI_Comm_size(PETSC_COMM_WORLD, &comm_size);
  if (comm_size != 1) {
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "affine_tet_stress_gram is a serial diagnostic; run with -np 1");
  }
 
  DMType dm_name = "DMMOFEM";
  CHKERR DMRegister_MoFEM(dm_name);
 
  moab::Core moab_core;
  moab::Interface &moab = moab_core;
  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
  if (pcomm == nullptr) {
    pcomm = new ParallelComm(&moab, PETSC_COMM_SELF);
  }
 
  EntityHandle meshset;
  CHKERR createAffineTet(moab, meshset);
 
  MoFEM::Core mofem_core(moab, PETSC_COMM_SELF, -2);
  MoFEM::Interface &m_field = mofem_core;
 
  BitRefLevel bit_level0 = BitRefLevel().set(0);
  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
      meshset, 3, bit_level0);
 
  EshelbianCore::hasNonHomogeneousMaterialBlock = false;
  EshelbianCore ep(m_field);
 
  ep.bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
  ep.bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
  ep.bcSpatialTractionVecPtr = boost::make_shared<TractionBcVec>();
  ep.bcSpatialAnalyticalTractionVecPtr =
      boost::make_shared<AnalyticalTractionBcVec>();
  ep.bcSpatialPressureVecPtr = boost::make_shared<PressureBcVec>();
  ep.bcSpatialNormalDisplacementVecPtr =
      boost::make_shared<NormalDisplacementBcVec>();
  ep.bcSpatialAnalyticalDisplacementVecPtr =
      boost::make_shared<AnalyticalDisplacementBcVec>();
  ep.externalStrainVecPtr = boost::make_shared<ExternalStrainVec>();
  ep.bcSpatialFreeTractionVecPtr = boost::make_shared<TractionFreeBc>();
  CHKERR ep.getTractionFreeBc(meshset, ep.bcSpatialFreeTractionVecPtr,
                              "CONTACT");
 
  CHKERR ep.addFields(meshset, include_bubble);
  CHKERR ep.projectGeometry(meshset, 0.0);
  CHKERR ep.addVolumeFiniteElement(meshset, include_bubble);
  CHKERR m_field.build_adjacencies(bit_level0, QUIET);
 
  auto dm = createDM(m_field.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "AFFINE_TET_EP", bit_level0,
                            BitRefLevel().set());
  CHKERR DMMoFEMSetDestroyProblem(dm, PETSC_TRUE);
  CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_FALSE);
  CHKERR DMMoFEMAddElement(dm, ep.elementVolumeName);
  m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
  CHKERR DMSetUp(dm);
  m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
 
  ep.solTSStep = createDMVector(dm);
  CHKERR VecZeroEntries(ep.solTSStep);
 
  auto diag_dm = createDM(m_field.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(diag_dm, dm, "AFFINE_TET_STRESS_DIAG");
  CHKERR DMMoFEMSetDestroyProblem(diag_dm, PETSC_TRUE);
  CHKERR DMMoFEMSetSquareProblem(diag_dm, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(diag_dm, PETSC_FALSE);
  CHKERR DMMoFEMAddElement(diag_dm, ep.elementVolumeName);
 
  switch (matrix_type) {
  case DiagnosticMatrixType::GRAM:
    CHKERR DMMoFEMAddSubFieldRow(diag_dm, ep.piolaStress);
    CHKERR DMMoFEMAddSubFieldCol(diag_dm, ep.piolaStress);
    if (include_bubble) {
      CHKERR DMMoFEMAddSubFieldRow(diag_dm, ep.bubbleField);
      CHKERR DMMoFEMAddSubFieldCol(diag_dm, ep.bubbleField);
    }
    break;
  case DiagnosticMatrixType::A00:
    if (!include_bubble) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "P-only diagnostic is implemented for -diag_matrix_type gram");
    }
    for (const auto field_name : {ep.stretchTensor, ep.bubbleField,
                                  ep.piolaStress, ep.rotAxis,
                                  ep.spatialL2Disp}) {
      CHKERR DMMoFEMAddSubFieldRow(diag_dm, field_name);
      CHKERR DMMoFEMAddSubFieldCol(diag_dm, field_name);
    }
    break;
  }
  CHKERR DMSetUp(diag_dm);
 
  CHKERR ep.addMaterial_Hencky(5.0, 0.25);
 
  auto fe_lhs = boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
  CHKERR ep.setBaseVolumeElementOps(adolc_tag, false, true, false, fe_lhs,
                                    include_bubble);
  switch (matrix_type) {
  case DiagnosticMatrixType::GRAM:
    if (include_bubble) {
      CHKERR ep.pushStressGramOps(fe_lhs);
    } else {
      CHKERR ep.pushPiolaStressGramOps(fe_lhs);
    }
    break;
  case DiagnosticMatrixType::A00:
    if (EshelbianCore::stretchHandling == EshelbianCore::NO_STREACH) {
      CHKERR ep.pushNoStretchVolumeA00Ops(fe_lhs);
    } else {
      CHKERR ep.pushStretchVolumeA00Ops(fe_lhs);
    }
    break;
  }
 
  auto A = createDMMatrix(diag_dm);
  CHKERR MatZeroEntries(A);
  fe_lhs->ksp_B = A;
  fe_lhs->ts_A = A;
  fe_lhs->ts_B = A;
  fe_lhs->ts_u = ep.solTSStep;
  fe_lhs->ts_u_t = ep.solTSStep;
  fe_lhs->ts_a = 0.0;
  fe_lhs->data_ctx = PetscData::CtxSetA | PetscData::CtxSetB |
                     PetscData::CtxSetX | PetscData::CtxSetX_T;
  CHKERR DMoFEMLoopFiniteElements(diag_dm, ep.elementVolumeName, fe_lhs);
  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 
  MatrixDouble dense;
  CHKERR getDenseMatrix(A, dense);
 
  CHKERR analyseDenseStressBlock(dense, matrix_type, include_bubble,
                                 rel_rank_tol, check_full_rank, view_matrix);
 
  MoFEMFunctionReturn(0);
}
 
/// Parses options, sets logging, and runs the diagnostic.
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  auto core_log = logging::core::get();
  core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "EP"));
  LogManager::setLog("EP");
  MOFEM_LOG_TAG("EP", "affine_tet_stress_gram");
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSelf(), "EPSELF"));
  LogManager::setLog("EPSELF");
  MOFEM_LOG_TAG("EPSELF", "affine_tet_stress_gram");
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSync(), "EPSYNC"));
  LogManager::setLog("EPSYNC");
  MOFEM_LOG_TAG("EPSYNC", "affine_tet_stress_gram");
 
  try {
    PetscBool view_matrix = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-diag_view_matrix",
                               &view_matrix, PETSC_NULLPTR);
    double rel_rank_tol = 1e-10;
    CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-diag_rank_tol",
                               &rel_rank_tol, PETSC_NULLPTR);
    PetscBool check_full_rank = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-diag_check_full_rank",
                               &check_full_rank, PETSC_NULLPTR);
    PetscBool include_bubble = PETSC_TRUE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-diag_include_bubble",
                               &include_bubble, PETSC_NULLPTR);
    DiagnosticMatrixType matrix_type;
    CHKERR getDiagnosticMatrixType(matrix_type);
    CHKERR assembleStressDiagnostic(matrix_type, include_bubble, rel_rank_tol,
                                    check_full_rank, view_matrix);
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
  return 0;
}