v0.8.23
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 1205 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 1206 of file EshelbianPlasticity.hpp.

1209  : OpAssembleVolume(row_field, data_ptr, OPROW), eps(eps) {
1210  sYmm = false;
1211 
1212 
1213 
1214  }
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 2101 of file EshelbianOperators.cpp.

2102  {
2104  if (row_type != MBTET)
2106 
2107  if (auto ep_fe_vol_ptr =
2108  dynamic_cast<const EpElement<VolumeElementForcesAndSourcesCore> *>(
2109  getFEMethod())) {
2110 
2111  AO aoSuu = PETSC_NULL;
2112  AO aoSBB = PETSC_NULL;
2113  AO aoSOO = PETSC_NULL;
2114  AO aoSww = PETSC_NULL;
2115 
2116  Mat Suu = PETSC_NULL;
2117  Mat SBB = PETSC_NULL;
2118  Mat SOO = PETSC_NULL;
2119  Mat Sww = PETSC_NULL;
2120 
2121  auto set_data = [&](auto fe) {
2122  aoSuu = fe->aoSuu;
2123  aoSBB = fe->aoSBubble;
2124  aoSOO = fe->aoSOmega;
2125  aoSww = fe->aoSw;
2126 
2127  Suu = fe->Suu;
2128  SBB = fe->SBubble;
2129  SOO = fe->SOmega;
2130  Sww = fe->Sw;
2131  };
2132  set_data(ep_fe_vol_ptr);
2133 
2134  if (Suu) {
2135 
2136  auto find_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2137  auto &row_name, auto &col_name, auto row_side,
2138  auto col_side, auto row_type, auto col_type) {
2139  return add_bmc.get<0>().find(boost::make_tuple(
2140  row_name, col_name, row_type, col_type, row_side, col_side));
2141  };
2142 
2143  auto set_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2144  auto &row_name, auto &col_name, auto row_side,
2145  auto col_side, auto row_type, auto col_type,
2146  const auto &m, const auto &row_ind,
2147  const auto &col_ind) {
2148  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2149  row_type, col_type);
2150  if (it != add_bmc.get<0>().end()) {
2151  it->setMat(m);
2152  it->setInd(row_ind, col_ind);
2153  it->setSetAtElement();
2154  return it;
2155  } else {
2156  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2157  row_name, col_name, row_type, col_type, row_side, col_side, m,
2158  row_ind, col_ind));
2159  return p_it.first;
2160  }
2161  };
2162 
2163  auto add_block = [&](DataAtIntegrationPts::BlockMatContainor &add_bmc,
2164  auto &row_name, auto &col_name, auto row_side,
2165  auto col_side, auto row_type, auto col_type,
2166  const auto &m, const auto &row_ind,
2167  const auto &col_ind) {
2168  auto it = find_block(add_bmc, row_name, col_name, row_side, col_side,
2169  row_type, col_type);
2170  if (it != add_bmc.get<0>().end()) {
2171  it->addMat(m);
2172  it->setInd(row_ind, col_ind);
2173  it->setSetAtElement();
2174  return it;
2175  } else {
2176  auto p_it = add_bmc.insert(DataAtIntegrationPts::BlockMatData(
2177  row_name, col_name, row_type, col_type, row_side, col_side, m,
2178  row_ind, col_ind));
2179  return p_it.first;
2180  }
2181  };
2182 
2183  auto assemble_block = [&](auto &bit, Mat S) {
2185  const VectorInt &rind = bit.rowInd;
2186  const VectorInt &cind = bit.colInd;
2187  const MatrixDouble &m = bit.M;
2188  CHKERR MatSetValues(S, rind.size(), &*rind.begin(), cind.size(),
2189  &*cind.begin(), &*m.data().begin(), ADD_VALUES);
2190 
2192  };
2193 
2194  auto invert_symm_mat = [&](MatrixDouble &m, auto &inv) {
2196  const int nb = m.size1();
2197 
2198  inv.resize(nb, nb, false);
2199  inv.clear();
2200  for (int cc = 0; cc != nb; ++cc)
2201  inv(cc, cc) = -1;
2202 
2203  iPIV.resize(nb, false);
2204  lapackWork.resize(nb * nb, false);
2205  const auto info = lapack_dsysv('L', nb, nb, &*m.data().begin(), nb,
2206  &*iPIV.begin(), &*inv.data().begin(), nb,
2207  &*lapackWork.begin(), nb * nb);
2208  if (info != 0)
2209  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2210  "Can not invert matrix info = %d", info);
2211 
2213  };
2214 
2215  auto invert_symm_schur = [&](DataAtIntegrationPts::BlockMatContainor &bmc,
2216  std::string field, auto &inv) {
2218 
2219  auto bit =
2220  bmc.get<1>().find(boost::make_tuple(field, field, MBTET, MBTET));
2221  if (bit != bmc.get<1>().end()) {
2222 
2223  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2224  CHKERR invert_symm_mat(m, inv);
2225 
2226  } else
2227  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2228  "%s matrix not found", field.c_str());
2229 
2231  };
2232 
2233  auto invert_nonsymm_mat = [&](MatrixDouble &m, auto &inv) {
2235 
2236  const int nb = m.size1();
2237 
2238  MatrixDouble trans_m = trans(m);
2239  MatrixDouble trans_inv;
2240  trans_inv.resize(nb, nb, false);
2241  trans_inv.clear();
2242  for (int c = 0; c != nb; ++c)
2243  trans_inv(c, c) = -1;
2244 
2245  iPIV.resize(nb, false);
2246  const auto info =
2247  lapack_dgesv(nb, nb, &*trans_m.data().begin(), nb, &*iPIV.begin(),
2248  &*trans_inv.data().begin(), nb);
2249  if (info != 0)
2250  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2251  "Can not invert matrix info = %d", info);
2252 
2253  inv.resize(nb, nb, false);
2254  noalias(inv) = trans(trans_inv);
2255 
2257  };
2258 
2259  auto invert_nonsymm_schur =
2260  [&](DataAtIntegrationPts::BlockMatContainor &bmc, std::string field,
2261  auto &inv, const bool debug = false) {
2263 
2264  auto bit = bmc.get<1>().find(
2265  boost::make_tuple(field, field, MBTET, MBTET));
2266  if (bit != bmc.get<1>().end()) {
2267 
2268  auto &m = *const_cast<MatrixDouble *>(&(bit->M));
2269  CHKERR invert_nonsymm_mat(m, inv);
2270 
2271  if (debug) {
2272  std::cerr << prod(m, inv) << endl;
2273  std::cerr << endl;
2274  }
2275 
2276  } else
2277  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2278  "%s matrix not found", field.c_str());
2279 
2281  };
2282 
2283  auto create_block_schur =
2286  std::string field, AO ao, MatrixDouble &inv) {
2288 
2289  for (auto &bit : add_bmc) {
2290  bit.unSetAtElement();
2291  bit.clearMat();
2292  }
2293 
2294  for (auto &bit : bmc) {
2295  if (bit.setAtElement && bit.rowField != field &&
2296  bit.colField != field) {
2297  VectorInt rind = bit.rowInd;
2298  VectorInt cind = bit.colInd;
2299  const MatrixDouble &m = bit.M;
2300  if (ao) {
2301  CHKERR AOApplicationToPetsc(ao, rind.size(), &*rind.begin());
2302  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2303  }
2304  auto it = set_block(add_bmc, bit.rowField, bit.colField,
2305  bit.rowSide, bit.colSide, bit.rowType,
2306  bit.colType, m, rind, cind);
2307  }
2308  }
2309 
2310  for (auto &bit_col : bmc) {
2311  if (bit_col.setAtElement && bit_col.rowField == field &&
2312  bit_col.colField != field) {
2313  const MatrixDouble &cm = bit_col.M;
2314  VectorInt cind = bit_col.colInd;
2315  invMat.resize(inv.size1(), cm.size2(), false);
2316  noalias(invMat) = prod(inv, cm);
2317  if (ao)
2318  CHKERR AOApplicationToPetsc(ao, cind.size(), &*cind.begin());
2319  for (auto &bit_row : bmc) {
2320  if (bit_row.setAtElement && bit_row.rowField != field &&
2321  bit_row.colField == field) {
2322  const MatrixDouble &rm = bit_row.M;
2323  VectorInt rind = bit_row.rowInd;
2324  K.resize(rind.size(), cind.size(), false);
2325  noalias(K) = prod(rm, invMat);
2326  if (ao)
2327  CHKERR AOApplicationToPetsc(ao, rind.size(),
2328  &*rind.begin());
2329  auto it = add_block(add_bmc, bit_row.rowField,
2330  bit_col.colField, bit_row.rowSide,
2331  bit_col.colSide, bit_row.rowType,
2332  bit_col.colType, K, rind, cind);
2333  }
2334  }
2335  }
2336  }
2337 
2339  };
2340 
2341  auto assemble_schur =
2342  [&](DataAtIntegrationPts::BlockMatContainor &add_bmc, Mat S,
2343  bool debug = false) {
2345  for (auto &bit : add_bmc) {
2346  if (bit.setAtElement)
2347  CHKERR assemble_block(bit, S);
2348  }
2349  if (debug) {
2350  for (auto &bit : add_bmc) {
2351  if (bit.setAtElement) {
2352  std::cerr << "assemble: " << bit.rowField << " "
2353  << bit.colField << endl;
2354  std::cerr << bit.M << endl;
2355  }
2356  }
2357  std::cerr << std::endl;
2358  }
2360  };
2361 
2362  auto precondition_schur =
2365  const std::string field, const MatrixDouble &diag_mat,
2366  const double eps) {
2368 
2369  for (auto &bit : add_bmc) {
2370  bit.unSetAtElement();
2371  bit.clearMat();
2372  }
2373 
2374  for (auto &bit : bmc) {
2375  if (bit.setAtElement) {
2376  if (bit.rowField != field || bit.colField != field)
2377  auto it =
2378  set_block(add_bmc, bit.rowField, bit.colField,
2379  bit.rowSide, bit.colSide, bit.rowType,
2380  bit.colType, bit.M, bit.rowInd, bit.colInd);
2381  }
2382  }
2383 
2384  auto bit = bmc.get<1>().find(
2385  boost::make_tuple(field, field, MBTET, MBTET));
2386  if (bit->setAtElement && bit != bmc.get<1>().end()) {
2387  auto it =
2388  set_block(add_bmc, bit->rowField, bit->colField, bit->rowSide,
2389  bit->colSide, bit->rowType, bit->colType, bit->M,
2390  bit->rowInd, bit->colInd);
2391  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2392  m += eps * diag_mat;
2393  } else {
2394  auto row_it = bmc.get<3>().lower_bound(field);
2395  for (; row_it != bmc.get<3>().upper_bound(field); ++row_it) {
2396  if (row_it->setAtElement) {
2397  auto it = set_block(add_bmc, field, field, 0, 0, MBTET, MBTET,
2398  diag_mat, row_it->rowInd, row_it->rowInd);
2399  MatrixDouble &m = const_cast<MatrixDouble &>(it->M);
2400  m *= eps;
2401  break;
2402  }
2403  }
2404  if (row_it == bmc.get<3>().end())
2405  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
2406  "row field not found %s", field.c_str());
2407  }
2408 
2410  };
2411 
2412  CHKERR invert_symm_schur(dataAtPts->blockMatContainor, "u",
2413  invBlockMat["uu"]);
2414  CHKERR create_block_schur(dataAtPts->blockMatContainor, blockMat["uu"],
2415  "u", aoSuu, invBlockMat["uu"]);
2416  CHKERR assemble_schur(blockMat["uu"], Suu);
2417 
2418  if (SBB) {
2419  CHKERR invert_symm_schur(blockMat["uu"], "BUBBLE", invBlockMat["BB"]);
2420  CHKERR create_block_schur(blockMat["uu"], blockMat["BB"], "BUBBLE",
2421  aoSBB, invBlockMat["BB"]);
2422  CHKERR precondition_schur(blockMat["BB"], blockMat["precBBOO"], "omega",
2423  *dataAtPts->ooMatPtr, eps);
2424  CHKERR assemble_schur(blockMat["precBBOO"], SBB);
2425 
2426  if (SOO) {
2427  CHKERR invert_nonsymm_schur(blockMat["precBBOO"], "omega",
2428  invBlockMat["OO"]);
2429  CHKERR create_block_schur(blockMat["precBBOO"], blockMat["OO"],
2430  "omega", aoSOO, invBlockMat["OO"]);
2431  CHKERR precondition_schur(blockMat["OO"], blockMat["precOOww"], "w",
2432  *dataAtPts->wwMatPtr, eps);
2433  CHKERR assemble_schur(blockMat["precOOww"], SOO);
2434 
2435  if (Sww) {
2436  CHKERR invert_symm_schur(blockMat["precOOww"], "w",
2437  invBlockMat["ww"]);
2438  CHKERR create_block_schur(blockMat["precOOww"], blockMat["ww"], "w",
2439  aoSww, invBlockMat["ww"]);
2440  CHKERR assemble_schur(blockMat["ww"], Sww);
2441  }
2442  }
2443  }
2444  }
2445  }
2446 
2448 } // namespace EshelbianPlasticity
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:500
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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:507
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:595
map< std::string, MatrixDouble > invBlockMat
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:75
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
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 1219 of file EshelbianPlasticity.hpp.

1219  {
1220  return assemble(side, type, data);
1221  }
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 1217 of file EshelbianPlasticity.hpp.

1217 { return 0; }

Member Data Documentation

◆ blockMat

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

Definition at line 1231 of file EshelbianPlasticity.hpp.

◆ eps

const double EshelbianPlasticity::OpSpatialSchurEnd::eps
private

Definition at line 1224 of file EshelbianPlasticity.hpp.

◆ invBlockMat

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

Definition at line 1230 of file EshelbianPlasticity.hpp.

◆ invMat

MatrixDouble EshelbianPlasticity::OpSpatialSchurEnd::invMat
private

Definition at line 1226 of file EshelbianPlasticity.hpp.

◆ iPIV

VectorInt EshelbianPlasticity::OpSpatialSchurEnd::iPIV
private

Definition at line 1227 of file EshelbianPlasticity.hpp.

◆ lapackWork

VectorDouble EshelbianPlasticity::OpSpatialSchurEnd::lapackWork
private

Definition at line 1228 of file EshelbianPlasticity.hpp.


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