v0.9.1
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 1137 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 1138 of file EshelbianPlasticity.hpp.

1141  : OpAssembleVolume(row_field, data_ptr, OPROW), eps(eps) {
1142  sYmm = false;
1143  }
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 2015 of file EshelbianOperators.cpp.

2016  {
2018  if (row_type != MBTET)
2020 
2021  if (auto ep_fe_vol_ptr =
2022  dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
2023  getFEMethod())) {
2024 
2025  AO aoSuu = PETSC_NULL;
2026  AO aoSBB = PETSC_NULL;
2027  AO aoSOO = PETSC_NULL;
2028  AO aoSww = PETSC_NULL;
2029 
2030  Mat Suu = PETSC_NULL;
2031  Mat SBB = PETSC_NULL;
2032  Mat SOO = PETSC_NULL;
2033  Mat Sww = PETSC_NULL;
2034 
2035  auto set_data = [&](auto fe) {
2036  aoSuu = fe->aoSuu;
2037  aoSBB = fe->aoSBubble;
2038  aoSOO = fe->aoSOmega;
2039  aoSww = fe->aoSw;
2040 
2041  Suu = fe->Suu;
2042  SBB = fe->SBubble;
2043  SOO = fe->SOmega;
2044  Sww = fe->Sw;
2045  };
2046  set_data(ep_fe_vol_ptr);
2047 
2048  if (Suu) {
2049 
2050  auto find_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2051  auto &row_name, auto &col_name, auto row_side,
2052  auto col_side, auto row_type, auto col_type) {
2053  return add_bmc.get<0>().find(boost::make_tuple(
2054  row_name, col_name, row_type, col_type, row_side, col_side));
2055  };
2056 
2057  auto set_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2058  auto &row_name, auto &col_name, auto row_side,
2059  auto col_side, auto row_type, auto col_type,
2060  const auto &m, const auto &row_ind,
2061  const auto &col_ind) {
2062  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2063  row_type, col_type);
2064  if (it != add_bmc.get<0>().end()) {
2065  it->setMat(m);
2066  it->setInd(row_ind, col_ind);
2067  it->setSetAtElement();
2068  return it;
2069  } else {
2070  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2071  row_name, col_name, row_type, col_type, row_side, col_side, m,
2072  row_ind, col_ind));
2073  return p_it.first;
2074  }
2075  };
2076 
2077  auto add_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2078  auto &row_name, auto &col_name, auto row_side,
2079  auto col_side, auto row_type, auto col_type,
2080  const auto &m, const auto &row_ind,
2081  const auto &col_ind) {
2082  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2083  row_type, col_type);
2084  if (it != add_bmc.get<0>().end()) {
2085  it->addMat(m);
2086  it->setInd(row_ind, col_ind);
2087  it->setSetAtElement();
2088  return it;
2089  } else {
2090  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2091  row_name, col_name, row_type, col_type, row_side, col_side, m,
2092  row_ind, col_ind));
2093  return p_it.first;
2094  }
2095  };
2096 
2097  auto assemble_block = [&](auto &bit, Mat S) {
2099  const VectorInt &rind = bit.rowInd;
2100  const VectorInt &cind = bit.colInd;
2101  const MatrixDouble &m = bit.M;
2102  CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
2103  &*cind.begin(), &*m.data().begin(), ADD_VALUES);
2104 
2106  };
2107 
2108  auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
2110  const int nb = m.size1();
2111 
2112  inv.resize(nb, nb, false);
2113  inv.clear();
2114  for (int cc = 0; cc != nb; ++cc)
2115  inv(cc, cc) = -1;
2116 
2117  iPIV.resize(nb, false);
2118  lapackWork.resize(nb * nb, false);
2119  const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
2120  &*iPIV.begin(), &*inv.data().begin(), nb,
2121  &*lapackWork.begin(), nb * nb);
2122  // if (info != 0)
2123  // SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2124  // "Can not invert matrix info = %d", info);
2125 
2127  };
2128 
2129  auto invert_symm_schur = [&](DataAtIntegrationPts::BlockMatContainor &bmc,
2130  std::string field, auto &inv) {
2132 
2133  auto bit =
2134  bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
2135  if (bit != bmc.get<1>().end()) {
2136 
2137  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2138  CHKERR invert_symm_mat(m, inv);
2139 
2140  } else
2141  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2142  "%s matrix not found", field.c_str());
2143 
2145  };
2146 
2147  auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
2149 
2150  const int nb = m.size1();
2151 
2152  MatrixDouble trans_m = trans(m);
2153  MatrixDouble trans_inv;
2154  trans_inv.resize(nb, nb, false);
2155  trans_inv.clear();
2156  for (int c = 0; c != nb; ++c)
2157  trans_inv(c, c) = -1;
2158 
2159  iPIV.resize(nb, false);
2160  const auto info =
2161  lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb, &*iPIV.begin(),
2162  &*trans_inv.data().begin(), nb);
2163  if (info != 0)
2164  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2165  "Can not invert matrix info = %d", info);
2166 
2167  inv.resize(nb, nb, false);
2168  noalias(inv) = trans(trans_inv);
2169 
2171  };
2172 
2173  auto invert_nonsymm_schur =
2174  [&](DataAtIntegrationPts::BlockMatContainor &bmc, std::string field,
2175  auto &inv, const bool debug = false) {
2177 
2178  auto bit = bmc.get<1>().find(
2179  boost::make_tuple(field, field, MBTET, MBTET));
2180  if (bit != bmc.get<1>().end()) {
2181 
2182  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2183  CHKERR invert_nonsymm_mat(m, inv);
2184 
2185  if (debug) {
2186  std::cerr << prod(m, inv) << endl;
2187  std::cerr << endl;
2188  }
2189 
2190  } else
2191  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2192  "%s matrix not found", field.c_str());
2193 
2195  };
2196 
2197  auto create_block_schur =
2200  std::string field, AO ao, MatrixDouble &inv) {
2202 
2203  for (auto &bit : add_bmc) {
2204  bit.unSetAtElement();
2205  bit.clearMat();
2206  }
2207 
2208  for (auto &bit : bmc) {
2209  if (bit.setAtElement && bit.rowField != field &&
2210  bit.colField != field) {
2211  VectorInt rind = bit.rowInd;
2212  VectorInt cind = bit.colInd;
2213  const MatrixDouble &m = bit.M;
2214  if (ao) {
2215  CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
2216  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2217  }
2218  auto it = set_block(add_bmc, bit.rowField, bit.colField,
2219  bit.rowSide, bit.colSide, bit.rowType,
2220  bit.colType, m, rind, cind);
2221  }
2222  }
2223 
2224  for (auto &bit_col : bmc) {
2225  if (bit_col.setAtElement && bit_col.rowField == field &&
2226  bit_col.colField != field) {
2227  const MatrixDouble &cm = bit_col.M;
2228  VectorInt cind = bit_col.colInd;
2229  invMat.resize(inv.size1(), cm.size2(), false);
2230  noalias(invMat) = prod(inv, cm);
2231  if (ao)
2232  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2233  for (auto &bit_row : bmc) {
2234  if (bit_row.setAtElement && bit_row.rowField != field &&
2235  bit_row.colField == field) {
2236  const MatrixDouble &rm = bit_row.M;
2237  VectorInt rind = bit_row.rowInd;
2238  K.resize(rind.size(), cind.size(), false);
2239  noalias(K) = prod(rm, invMat);
2240  if (ao)
2241  CHKERR AOApplicationToPetsc(ao, rind.size(),
2242  &*rind.begin());
2243  auto it = add_block(add_bmc, bit_row.rowField,
2244  bit_col.colField, bit_row.rowSide,
2245  bit_col.colSide, bit_row.rowType,
2246  bit_col.colType, K, rind, cind);
2247  }
2248  }
2249  }
2250  }
2251 
2253  };
2254 
2255  auto assemble_schur =
2256  [&](DataAtIntegrationPts::BlockMatContainor &add_bmc, Mat S,
2257  bool debug = false) {
2259  for (auto &bit : add_bmc) {
2260  if (bit.setAtElement)
2261  CHKERR assemble_block(bit, S);
2262  }
2263  if (debug) {
2264  for (auto &bit : add_bmc) {
2265  if (bit.setAtElement) {
2266  std::cerr << "assemble: " << bit.rowField << " "
2267  << bit.colField << endl;
2268  std::cerr << bit.M << endl;
2269  }
2270  }
2271  std::cerr << std::endl;
2272  }
2274  };
2275 
2276  auto precondition_schur =
2279  const std::string field, const MatrixDouble &diag_mat,
2280  const double eps) {
2282 
2283  for (auto &bit : add_bmc) {
2284  bit.unSetAtElement();
2285  bit.clearMat();
2286  }
2287 
2288  for (auto &bit : bmc) {
2289  if (bit.setAtElement) {
2290  if (bit.rowField != field || bit.colField != field)
2291  auto it =
2292  set_block(add_bmc, bit.rowField, bit.colField,
2293  bit.rowSide, bit.colSide, bit.rowType,
2294  bit.colType, bit.M, bit.rowInd, bit.colInd);
2295  }
2296  }
2297 
2298  auto bit = bmc.get<1>().find(
2299  boost::make_tuple(field, field, MBTET, MBTET));
2300  if (bit->setAtElement && bit != bmc.get<1>().end()) {
2301  auto it =
2302  set_block(add_bmc, bit->rowField, bit->colField, bit->rowSide,
2303  bit->colSide, bit->rowType, bit->colType, bit->M,
2304  bit->rowInd, bit->colInd);
2305  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2306  m += eps * diag_mat;
2307  } else {
2308  auto row_it = bmc.get<3>().lower_bound(field);
2309  for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
2310  if (row_it->setAtElement) {
2311  auto it = set_block(add_bmc, field, field, 0, 0, MBTET, MBTET,
2312  diag_mat, row_it->rowInd, row_it->rowInd);
2313  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2314  m *= eps;
2315  break;
2316  }
2317  }
2318  if (row_it == bmc.get<3>().end())
2319  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2320  "row field not found %s", field.c_str());
2321  }
2322 
2324  };
2325 
2326  CHKERR invert_symm_schur(dataAtPts->blockMatContainor, "u",
2327  invBlockMat["uu"]);
2328  CHKERR create_block_schur(dataAtPts->blockMatContainor, blockMat["uu"],
2329  "u", aoSuu, invBlockMat["uu"]);
2330  CHKERR assemble_schur(blockMat["uu"], Suu);
2331 
2332  if (SBB) {
2333  CHKERR invert_symm_schur(blockMat["uu"], "BUBBLE", invBlockMat["BB"]);
2334  CHKERR create_block_schur(blockMat["uu"], blockMat["BB"], "BUBBLE",
2335  aoSBB, invBlockMat["BB"]);
2336  CHKERR precondition_schur(blockMat["BB"], blockMat["precBBOO"], "omega",
2337  *dataAtPts->ooMatPtr, eps);
2338  CHKERR assemble_schur(blockMat["precBBOO"], SBB);
2339 
2340  if (SOO) {
2341  CHKERR invert_nonsymm_schur(blockMat["precBBOO"], "omega",
2342  invBlockMat["OO"]);
2343  CHKERR create_block_schur(blockMat["precBBOO"], blockMat["OO"],
2344  "omega", aoSOO, invBlockMat["OO"]);
2345  if (dataAtPts->wwMatPtr) {
2346  CHKERR precondition_schur(blockMat["OO"], blockMat["precOOww"], "w",
2347  *dataAtPts->wwMatPtr, -eps);
2348  } else {
2349  blockMat["precOOww"] = blockMat["OO"];
2350  }
2351  CHKERR assemble_schur(blockMat["precOOww"], SOO);
2352 
2353  if (Sww) {
2354  CHKERR invert_symm_schur(blockMat["precOOww"], "w",
2355  invBlockMat["ww"]);
2356  CHKERR create_block_schur(blockMat["precOOww"], blockMat["ww"], "w",
2357  aoSww, invBlockMat["ww"]);
2358  CHKERR assemble_schur(blockMat["ww"], Sww);
2359  }
2360  }
2361  }
2362  }
2363  }
2364 
2366 }
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:506
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
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:513
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
static const bool debug
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
#define CHKERR
Inline error check.
Definition: definitions.h:601
map< std::string, MatrixDouble > invBlockMat
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:71
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
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 1147 of file EshelbianPlasticity.hpp.

1147  {
1148  return assemble(side, type, data);
1149  }
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 1145 of file EshelbianPlasticity.hpp.

1145 { return 0; }

Member Data Documentation

◆ blockMat

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

Definition at line 1159 of file EshelbianPlasticity.hpp.

◆ eps

const double EshelbianPlasticity::OpSpatialSchurEnd::eps
private

Definition at line 1152 of file EshelbianPlasticity.hpp.

◆ invBlockMat

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

Definition at line 1158 of file EshelbianPlasticity.hpp.

◆ invMat

MatrixDouble EshelbianPlasticity::OpSpatialSchurEnd::invMat
private

Definition at line 1154 of file EshelbianPlasticity.hpp.

◆ iPIV

VectorInt EshelbianPlasticity::OpSpatialSchurEnd::iPIV
private

Definition at line 1155 of file EshelbianPlasticity.hpp.

◆ lapackWork

VectorDouble EshelbianPlasticity::OpSpatialSchurEnd::lapackWork
private

Definition at line 1156 of file EshelbianPlasticity.hpp.


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