v0.10.0
Public Member Functions | Private Attributes | List of all members
EshelbianPlasticity::OpSpatialSchurEnd Struct Reference

#include <users_modules/eshelbian_plasticty/src/EshelbianPlasticity.hpp>

Inheritance diagram for EshelbianPlasticity::OpSpatialSchurEnd:
[legend]
Collaboration diagram for EshelbianPlasticity::OpSpatialSchurEnd:
[legend]

Public Member Functions

 OpSpatialSchurEnd (const std::string &row_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const double eps)
 
MoFEMErrorCode integrate (EntData &row_data)
 
MoFEMErrorCode assemble (int side, EntityType type, EntData &data)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 
- Public Member Functions inherited from EshelbianPlasticity::OpAssembleVolume
 OpAssembleVolume (const std::string &field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type)
 
 OpAssembleVolume (const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type, const bool assemble_symmetry)
 
- Public Member Functions inherited from EshelbianPlasticity::OpAssembleBasic< VolUserDataOperator >
 OpAssembleBasic (const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type)
 
 OpAssembleBasic (const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type, const bool assemble_symmetry)
 
virtual MoFEMErrorCode integrate (int row_side, EntityType row_type, EntData &data)
 
virtual MoFEMErrorCode integrate (EntData &row_data, EntData &col_data)
 
virtual MoFEMErrorCode assemble (EntData &data)
 
virtual MoFEMErrorCode assemble (int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 
MoFEMErrorCode doWork (int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
 

Private Attributes

const double eps
 
MatrixDouble invMat
 
VectorInt iPIV
 
VectorDouble lapackWork
 
map< std::string, MatrixDouble > invBlockMat
 
map< std::string, DataAtIntegrationPts::BlockMatContainorblockMat
 

Additional Inherited Members

- Public Attributes inherited from EshelbianPlasticity::OpAssembleBasic< VolUserDataOperator >
const bool assembleSymmetry
 
boost::shared_ptr< DataAtIntegrationPtsdataAtPts
 data at integration pts More...
 
VectorDouble nF
 local right hand side vector More...
 
MatrixDouble K
 local tangent matrix More...
 
MatrixDouble transposeK
 

Detailed Description

Definition at line 1145 of file EshelbianPlasticity.hpp.

Constructor & Destructor Documentation

◆ OpSpatialSchurEnd()

EshelbianPlasticity::OpSpatialSchurEnd::OpSpatialSchurEnd ( const std::string &  row_field,
boost::shared_ptr< DataAtIntegrationPts > &  data_ptr,
const double  eps 
)

Definition at line 1146 of file EshelbianPlasticity.hpp.

1149  : OpAssembleVolume(row_field, data_ptr, OPROW), eps(eps) {
1150  sYmm = false;
1151  }
OpAssembleVolume(const std::string &field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type)

Member Function Documentation

◆ assemble()

MoFEMErrorCode EshelbianPlasticity::OpSpatialSchurEnd::assemble ( int  side,
EntityType  type,
EntData data 
)
virtual

Reimplemented from EshelbianPlasticity::OpAssembleBasic< VolUserDataOperator >.

Definition at line 1950 of file EshelbianOperators.cpp.

1951  {
1953  if (row_type != MBTET)
1955 
1956  if (auto ep_fe_vol_ptr =
1957  dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
1958  getFEMethod())) {
1959 
1960  AO aoSuu = PETSC_NULL;
1961  AO aoSBB = PETSC_NULL;
1962  AO aoSOO = PETSC_NULL;
1963  AO aoSww = PETSC_NULL;
1964 
1965  Mat Suu = PETSC_NULL;
1966  Mat SBB = PETSC_NULL;
1967  Mat SOO = PETSC_NULL;
1968  Mat Sww = PETSC_NULL;
1969 
1970  auto set_data = [&](auto fe) {
1971  aoSuu = fe->aoSuu;
1972  aoSBB = fe->aoSBubble;
1973  aoSOO = fe->aoSOmega;
1974  aoSww = fe->aoSw;
1975 
1976  Suu = fe->Suu;
1977  SBB = fe->SBubble;
1978  SOO = fe->SOmega;
1979  Sww = fe->Sw;
1980  };
1981  set_data(ep_fe_vol_ptr);
1982 
1983  if (Suu) {
1984 
1985  auto find_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1986  auto &row_name, auto &col_name, auto row_side,
1987  auto col_side, auto row_type, auto col_type) {
1988  return add_bmc.get<0>().find(boost::make_tuple(
1989  row_name, col_name, row_type, col_type, row_side, col_side));
1990  };
1991 
1992  auto set_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
1993  auto &row_name, auto &col_name, auto row_side,
1994  auto col_side, auto row_type, auto col_type,
1995  const auto &m, const auto &row_ind,
1996  const auto &col_ind) {
1997  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
1998  row_type, col_type);
1999  if (it != add_bmc.get<0>().end()) {
2000  it->setMat(m);
2001  it->setInd(row_ind, col_ind);
2002  it->setSetAtElement();
2003  return it;
2004  } else {
2005  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2006  row_name, col_name, row_type, col_type, row_side, col_side, m,
2007  row_ind, col_ind));
2008  return p_it.first;
2009  }
2010  };
2011 
2012  auto add_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2013  auto &row_name, auto &col_name, auto row_side,
2014  auto col_side, auto row_type, auto col_type,
2015  const auto &m, const auto &row_ind,
2016  const auto &col_ind) {
2017  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2018  row_type, col_type);
2019  if (it != add_bmc.get<0>().end()) {
2020  it->addMat(m);
2021  it->setInd(row_ind, col_ind);
2022  it->setSetAtElement();
2023  return it;
2024  } else {
2025  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2026  row_name, col_name, row_type, col_type, row_side, col_side, m,
2027  row_ind, col_ind));
2028  return p_it.first;
2029  }
2030  };
2031 
2032  auto assemble_block = [&](auto &bit, Mat S) {
2034  const VectorInt &rind = bit.rowInd;
2035  const VectorInt &cind = bit.colInd;
2036  const MatrixDouble &m = bit.M;
2037  CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
2038  &*cind.begin(), &*m.data().begin(), ADD_VALUES);
2039 
2041  };
2042 
2043  auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
2045  const int nb = m.size1();
2046 
2047  inv.resize(nb, nb, false);
2048  inv.clear();
2049  for (int cc = 0; cc != nb; ++cc)
2050  inv(cc, cc) = -1;
2051 
2052  iPIV.resize(nb, false);
2053  lapackWork.resize(nb * nb, false);
2054  const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
2055  &*iPIV.begin(), &*inv.data().begin(), nb,
2056  &*lapackWork.begin(), nb * nb);
2057  // if (info != 0)
2058  // SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2059  // "Can not invert matrix info = %d", info);
2060 
2062  };
2063 
2064  auto invert_symm_schur = [&](DataAtIntegrationPts::BlockMatContainor &bmc,
2065  std::string field, auto &inv) {
2067 
2068  auto bit =
2069  bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
2070  if (bit != bmc.get<1>().end()) {
2071 
2072  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2073  CHKERR invert_symm_mat(m, inv);
2074 
2075  } else
2076  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2077  "%s matrix not found", field.c_str());
2078 
2080  };
2081 
2082  auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
2084 
2085  const int nb = m.size1();
2086 
2087  MatrixDouble trans_m = trans(m);
2088  MatrixDouble trans_inv;
2089  trans_inv.resize(nb, nb, false);
2090  trans_inv.clear();
2091  for (int c = 0; c != nb; ++c)
2092  trans_inv(c, c) = -1;
2093 
2094  iPIV.resize(nb, false);
2095  const auto info =
2096  lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb, &*iPIV.begin(),
2097  &*trans_inv.data().begin(), nb);
2098  if (info != 0)
2099  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2100  "Can not invert matrix info = %d", info);
2101 
2102  inv.resize(nb, nb, false);
2103  noalias(inv) = trans(trans_inv);
2104 
2106  };
2107 
2108  auto invert_nonsymm_schur =
2109  [&](DataAtIntegrationPts::BlockMatContainor &bmc, std::string field,
2110  auto &inv, const bool debug = false) {
2112 
2113  auto bit = bmc.get<1>().find(
2114  boost::make_tuple(field, field, MBTET, MBTET));
2115  if (bit != bmc.get<1>().end()) {
2116 
2117  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2118  CHKERR invert_nonsymm_mat(m, inv);
2119 
2120  if (debug) {
2121  std::cerr << prod(m, inv) << endl;
2122  std::cerr << endl;
2123  }
2124 
2125  } else
2126  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2127  "%s matrix not found", field.c_str());
2128 
2130  };
2131 
2132  auto create_block_schur =
2135  std::string field, AO ao, MatrixDouble &inv) {
2137 
2138  for (auto &bit : add_bmc) {
2139  bit.unSetAtElement();
2140  bit.clearMat();
2141  }
2142 
2143  for (auto &bit : bmc) {
2144  if (bit.setAtElement && bit.rowField != field &&
2145  bit.colField != field) {
2146  VectorInt rind = bit.rowInd;
2147  VectorInt cind = bit.colInd;
2148  const MatrixDouble &m = bit.M;
2149  if (ao) {
2150  CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
2151  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2152  }
2153  auto it = set_block(add_bmc, bit.rowField, bit.colField,
2154  bit.rowSide, bit.colSide, bit.rowType,
2155  bit.colType, m, rind, cind);
2156  }
2157  }
2158 
2159  for (auto &bit_col : bmc) {
2160  if (bit_col.setAtElement && bit_col.rowField == field &&
2161  bit_col.colField != field) {
2162  const MatrixDouble &cm = bit_col.M;
2163  VectorInt cind = bit_col.colInd;
2164  invMat.resize(inv.size1(), cm.size2(), false);
2165  noalias(invMat) = prod(inv, cm);
2166  if (ao)
2167  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2168  for (auto &bit_row : bmc) {
2169  if (bit_row.setAtElement && bit_row.rowField != field &&
2170  bit_row.colField == field) {
2171  const MatrixDouble &rm = bit_row.M;
2172  VectorInt rind = bit_row.rowInd;
2173  K.resize(rind.size(), cind.size(), false);
2174  noalias(K) = prod(rm, invMat);
2175  if (ao)
2176  CHKERR AOApplicationToPetsc(ao, rind.size(),
2177  &*rind.begin());
2178  auto it = add_block(add_bmc, bit_row.rowField,
2179  bit_col.colField, bit_row.rowSide,
2180  bit_col.colSide, bit_row.rowType,
2181  bit_col.colType, K, rind, cind);
2182  }
2183  }
2184  }
2185  }
2186 
2188  };
2189 
2190  auto assemble_schur =
2191  [&](DataAtIntegrationPts::BlockMatContainor &add_bmc, Mat S,
2192  bool debug = false) {
2194  for (auto &bit : add_bmc) {
2195  if (bit.setAtElement)
2196  CHKERR assemble_block(bit, S);
2197  }
2198  if (debug) {
2199  for (auto &bit : add_bmc) {
2200  if (bit.setAtElement) {
2201  std::cerr << "assemble: " << bit.rowField << " "
2202  << bit.colField << endl;
2203  std::cerr << bit.M << endl;
2204  }
2205  }
2206  std::cerr << std::endl;
2207  }
2209  };
2210 
2211  auto precondition_schur =
2214  const std::string field, const MatrixDouble &diag_mat,
2215  const double eps) {
2217 
2218  for (auto &bit : add_bmc) {
2219  bit.unSetAtElement();
2220  bit.clearMat();
2221  }
2222 
2223  for (auto &bit : bmc) {
2224  if (bit.setAtElement) {
2225  if (bit.rowField != field || bit.colField != field)
2226  auto it =
2227  set_block(add_bmc, bit.rowField, bit.colField,
2228  bit.rowSide, bit.colSide, bit.rowType,
2229  bit.colType, bit.M, bit.rowInd, bit.colInd);
2230  }
2231  }
2232 
2233  auto bit = bmc.get<1>().find(
2234  boost::make_tuple(field, field, MBTET, MBTET));
2235  if (bit->setAtElement && bit != bmc.get<1>().end()) {
2236  auto it =
2237  set_block(add_bmc, bit->rowField, bit->colField, bit->rowSide,
2238  bit->colSide, bit->rowType, bit->colType, bit->M,
2239  bit->rowInd, bit->colInd);
2240  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2241  m += eps * diag_mat;
2242  } else {
2243  auto row_it = bmc.get<3>().lower_bound(field);
2244  for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
2245  if (row_it->setAtElement) {
2246  auto it = set_block(add_bmc, field, field, 0, 0, MBTET, MBTET,
2247  diag_mat, row_it->rowInd, row_it->rowInd);
2248  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2249  m *= eps;
2250  break;
2251  }
2252  }
2253  if (row_it == bmc.get<3>().end())
2254  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2255  "row field not found %s", field.c_str());
2256  }
2257 
2259  };
2260 
2261  CHKERR invert_symm_schur(dataAtPts->blockMatContainor, "u",
2262  invBlockMat["uu"]);
2263  CHKERR create_block_schur(dataAtPts->blockMatContainor, blockMat["uu"],
2264  "u", aoSuu, invBlockMat["uu"]);
2265  CHKERR assemble_schur(blockMat["uu"], Suu);
2266 
2267  if (SBB) {
2268  CHKERR invert_symm_schur(blockMat["uu"], "BUBBLE", invBlockMat["BB"]);
2269  CHKERR create_block_schur(blockMat["uu"], blockMat["BB"], "BUBBLE",
2270  aoSBB, invBlockMat["BB"]);
2271  CHKERR precondition_schur(blockMat["BB"], blockMat["precBBOO"], "omega",
2272  *dataAtPts->ooMatPtr, eps);
2273  CHKERR assemble_schur(blockMat["precBBOO"], SBB);
2274 
2275  if (SOO) {
2276  CHKERR invert_nonsymm_schur(blockMat["precBBOO"], "omega",
2277  invBlockMat["OO"]);
2278  CHKERR create_block_schur(blockMat["precBBOO"], blockMat["OO"],
2279  "omega", aoSOO, invBlockMat["OO"]);
2280  if (dataAtPts->wwMatPtr) {
2281  CHKERR precondition_schur(blockMat["OO"], blockMat["precOOww"], "w",
2282  *dataAtPts->wwMatPtr, -eps);
2283  } else {
2284  blockMat["precOOww"] = blockMat["OO"];
2285  }
2286  CHKERR assemble_schur(blockMat["precOOww"], SOO);
2287 
2288  if (Sww) {
2289  CHKERR invert_symm_schur(blockMat["precOOww"], "w",
2290  invBlockMat["ww"]);
2291  CHKERR create_block_schur(blockMat["precOOww"], blockMat["ww"], "w",
2292  aoSww, invBlockMat["ww"]);
2293  CHKERR assemble_schur(blockMat["ww"], Sww);
2294  }
2295  }
2296  }
2297  }
2298  }
2299 
2301 }
map< std::string, DataAtIntegrationPts::BlockMatContainor > blockMat
multi_index_container< BlockMatData, indexed_by< ordered_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType >, member< BlockMatData, int, &BlockMatData::rowSide >, member< BlockMatData, int, &BlockMatData::colSide > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField > > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::rowField > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::colField > > > > BlockMatContainor
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
static Index< 'm', 3 > m
static const bool debug
#define CHKERR
Inline error check.
Definition: definitions.h:604
map< std::string, MatrixDouble > invBlockMat
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
static __CLPK_integer lapack_dsysv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:232
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:188

◆ doWork()

MoFEMErrorCode EshelbianPlasticity::OpSpatialSchurEnd::doWork ( int  side,
EntityType  type,
EntData data 
)

Definition at line 1155 of file EshelbianPlasticity.hpp.

1155  {
1156  return assemble(side, type, data);
1157  }
MoFEMErrorCode assemble(int side, EntityType type, EntData &data)

◆ integrate()

MoFEMErrorCode EshelbianPlasticity::OpSpatialSchurEnd::integrate ( EntData row_data)
virtual

Reimplemented from EshelbianPlasticity::OpAssembleBasic< VolUserDataOperator >.

Definition at line 1153 of file EshelbianPlasticity.hpp.

1153 { return 0; }

Member Data Documentation

◆ blockMat

map<std::string, DataAtIntegrationPts::BlockMatContainor> EshelbianPlasticity::OpSpatialSchurEnd::blockMat
private

Definition at line 1167 of file EshelbianPlasticity.hpp.

◆ eps

const double EshelbianPlasticity::OpSpatialSchurEnd::eps
private

Definition at line 1160 of file EshelbianPlasticity.hpp.

◆ invBlockMat

map<std::string, MatrixDouble> EshelbianPlasticity::OpSpatialSchurEnd::invBlockMat
private

Definition at line 1166 of file EshelbianPlasticity.hpp.

◆ invMat

MatrixDouble EshelbianPlasticity::OpSpatialSchurEnd::invMat
private

Definition at line 1162 of file EshelbianPlasticity.hpp.

◆ iPIV

VectorInt EshelbianPlasticity::OpSpatialSchurEnd::iPIV
private

Definition at line 1163 of file EshelbianPlasticity.hpp.

◆ lapackWork

VectorDouble EshelbianPlasticity::OpSpatialSchurEnd::lapackWork
private

Definition at line 1164 of file EshelbianPlasticity.hpp.


The documentation for this struct was generated from the following files: