@@ -1670,20 +1670,22 @@ _computeK(const Real& lambda, const Real& mu, const DoFLocalId& igauss, const In
16701670 VariableDoFReal& gauss_jacobian (m_gauss_on_cells.gaussJacobian ());
16711671
16721672 auto jacobian = gauss_jacobian[igauss];
1673-
1673+ auto a{ lambda + 2 .*mu };
16741674 auto size{NDIM * nb_nodes};
16751675
16761676 // Setting the "B" matrix size for the max number of nodes in 3D:
16771677 // 8 nodes for a lin element/20 nodes for a quadratic one
1678- RealUniqueArray2 Bmat (NDIM, size );
1678+ RealUniqueArray2 Bmat = _getB (igauss,nb_nodes );
16791679
1680- auto a{ lambda + 2 .*mu };
1680+ /*
1681+ RealUniqueArray2 Bmat(NDIM, size);
16811682
16821683 for (int i = 0; i < NDIM; ++i)
16831684 for (int j = 0; j < size; ++j) {
16841685 Bmat(i, j) = 0.;
16851686 }
16861687
1688+
16871689 // ! Computes the Inverse Jacobian Matrix of a 2D or 3D finite-element
16881690 auto jac = gauss_jacobmat[igauss];
16891691 Real3x3 ijac;
@@ -1697,22 +1699,21 @@ _computeK(const Real& lambda, const Real& mu, const DoFLocalId& igauss, const In
16971699 ijac.y.x = -jac.y.x / jacobian;
16981700 ijac.y.y = jac.x.x / jacobian;
16991701 }
1700-
1701- auto wt = gauss_weight[igauss] * jacobian;
1702-
1703- // ------------------------------------------------------
1704- // Elementary Derivation Matrix B at current Gauss point
1705- // ------------------------------------------------------
1706- for (Int32 inod = 0 ; inod < nb_nodes; ++inod) {
1707- auto dPhi = gauss_shapederiv[igauss][inod];
1708- for (int i = 0 ; i < NDIM; ++i){
1709- auto bi{0 .};
1710- for (int j = 0 ; j < NDIM; ++j) {
1711- bi += ijac[i][j] * dPhi[j];
1712- }
1713- Bmat (i, inod) = bi;
1714- }
1702+ //------------------------------------------------------
1703+ // Elementary Derivation Matrix B at current Gauss point
1704+ //------------------------------------------------------
1705+ for (Int32 inod = 0; inod < nb_nodes; ++inod) {
1706+ auto dPhi = gauss_shapederiv[igauss][inod];
1707+ for (int i = 0; i < NDIM; ++i){
1708+ auto bi{0.};
1709+ for (int j = 0; j < NDIM; ++j) {
1710+ bi += ijac[i][j] * dPhi[j];
1711+ }
1712+ Bmat(i, inod) = bi;
17151713 }
1714+ }
1715+ */
1716+ auto wt = gauss_weight[igauss] * jacobian;
17161717
17171718 // ----------------------------------------------
17181719 // Elementary Stiffness (Ke) Matrix assembly
@@ -1803,7 +1804,12 @@ _computeK(const Real& lambda, const Real& mu, const DoFLocalId& igauss, const In
18031804 D(0,1) = D(0,2) = D(1,2) = lambda
18041805 All other terms = 0.
18051806------------------------------------------------------------------------------------*/
1806- auto kij = wt * (Bii (0 ) * (a * Bjj (0 ) + lambda * Bjj (1 ) + lambda * Bjj (2 )) + Bii (1 ) * (lambda * Bjj (0 ) + a * Bjj (1 ) + lambda * Bjj (2 )) + Bii (2 ) * (lambda * Bjj (0 ) + lambda * Bjj (1 ) + a * Bjj (2 )) + Bii (3 ) * (mu * Bjj (3 )) + Bii (4 ) * (mu * Bjj (4 )) + Bii (5 ) * (mu * Bjj (5 )));
1807+ auto kij = wt * (Bii (0 ) * (a * Bjj (0 ) + lambda * Bjj (1 ) + lambda * Bjj (2 ))
1808+ + Bii (1 ) * (lambda * Bjj (0 ) + a * Bjj (1 ) + lambda * Bjj (2 ))
1809+ + Bii (2 ) * (lambda * Bjj (0 ) + lambda * Bjj (1 ) + a * Bjj (2 ))
1810+ + Bii (3 ) * (mu * Bjj (3 ))
1811+ + Bii (4 ) * (mu * Bjj (4 ))
1812+ + Bii (5 ) * (mu * Bjj (5 )));
18071813
18081814 Ke (ii, jj) = kij;
18091815 Ke (jj, ii) = kij;
@@ -2340,7 +2346,60 @@ _assembleLinearRHS(){
23402346 }
23412347 }
23422348}
2349+ /* ---------------------------------------------------------------------------*/
2350+ /* ---------------------------------------------------------------------------*/
2351+ // ! Compute Elementary Derivation Matrix B at current Gauss point
2352+
2353+ RealUniqueArray2 NLDynamicModule::
2354+ _getB (const DoFLocalId& igauss, const Int32& nb_nodes)
2355+ {
23432356
2357+ auto gauss_point (m_gauss_on_cells.gaussCellConnectivityView ());
2358+ VariableDoFArrayReal3& gauss_shapederiv (m_gauss_on_cells.gaussShapeDeriv ());
2359+ VariableDoFReal3x3& gauss_jacobmat (m_gauss_on_cells.gaussJacobMat ());
2360+ VariableDoFReal& gauss_jacobian (m_gauss_on_cells.gaussJacobian ());
2361+
2362+ auto jacobian = gauss_jacobian[igauss];
2363+ auto size{ NDIM * nb_nodes };
2364+
2365+ // Setting the "B" matrix size for the max number of nodes in 3D:
2366+ // 8 nodes for a lin element/20 nodes for a quadratic one
2367+ RealUniqueArray2 Bmat (NDIM, size);
2368+
2369+ for (int i = 0 ; i < NDIM; ++i)
2370+ for (int j = 0 ; j < size; ++j) {
2371+ Bmat (i, j) = 0 .;
2372+ }
2373+
2374+ // ! Computes the Inverse Jacobian Matrix of a 2D or 3D finite-element
2375+ auto jac = gauss_jacobmat[igauss];
2376+ Real3x3 ijac;
2377+
2378+ if (NDIM == 3 ) {
2379+ ijac = math::inverseMatrix (jac);
2380+ }
2381+ else {
2382+ ijac.x .x = jac.y .y / jacobian;
2383+ ijac.x .y = -jac.x .y / jacobian;
2384+ ijac.y .x = -jac.y .x / jacobian;
2385+ ijac.y .y = jac.x .x / jacobian;
2386+ }
2387+
2388+ // ------------------------------------------------------
2389+ // Elementary Derivation Matrix B at current Gauss point
2390+ // ------------------------------------------------------
2391+ for (Int32 inod = 0 ; inod < nb_nodes; ++inod) {
2392+ auto dPhi = gauss_shapederiv[igauss][inod];
2393+ for (int i = 0 ; i < NDIM; ++i) {
2394+ auto bi{ 0 . };
2395+ for (int j = 0 ; j < NDIM; ++j) {
2396+ bi += ijac[i][j] * dPhi[j];
2397+ }
2398+ Bmat (i, inod) = bi;
2399+ }
2400+ }
2401+ return Bmat;
2402+ }
23442403/* ---------------------------------------------------------------------------*/
23452404/* ---------------------------------------------------------------------------*/
23462405// ! Stress prediction
@@ -2349,6 +2408,7 @@ void NLDynamicModule::stress_prediction(bool init, bool isRef)
23492408{
23502409 auto gauss_point (m_gauss_on_cells.gaussCellConnectivityView ());
23512410
2411+ VariableDoFReal& gauss_weight (m_gauss_on_cells.gaussWeight ());
23522412 VariableDoFArrayReal3x3& gauss_stress (m_gauss_on_cells.gaussStress ());
23532413 VariableDoFArrayReal3x3& gauss_strain (m_gauss_on_cells.gaussStrain ());
23542414 VariableDoFArrayReal3x3& gauss_strain_plastic (m_gauss_on_cells.gaussStrainPlastic ());
@@ -2378,9 +2438,10 @@ void NLDynamicModule::stress_prediction(bool init, bool isRef)
23782438 sign.fromReal3x3ToTensor2 (gauss_stress[gauss_pti][1 ]);
23792439 epspn.fromReal3x3ToTensor2 (gauss_strain_plastic[gauss_pti][1 ]);
23802440
2441+ // Compute strain increment for this interation
23812442 Tensor2 deps;
23822443 {
2383- // calcul deps à faire avec unodes
2444+ RealUniqueArray2 Bmat = _getB (gauss_pti,cell_nbnod);
23842445 }
23852446
23862447 for (Int32 ip = 0 ; ip < nblaw; ++ip) {
0 commit comments