2D struct simulation fails

[ktbolt]

Description

A simulation using the struct equation with the nHK material model fails for 2D triangle elements

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 ST 1-1  1.100e-01  [0 1.000e+00 1.000e+00 9.405e-05]  [296 -1 54]
 ST 1-2  2.670e-01  [-3 6.816e-01 6.816e-01 9.201e-05]  [572 -4 71]
 ST 1-3  5.140e-01  !93 4.841e+04 4.841e+04 2.385e-01!  !1000 0 82!  WARNING: The linear system solution has not converged
 ST 1-4  7.550e-01  !725 1.844e+36 1.844e+36 1.672e+02!  !1000 0 81!  WARNING: The linear system solution has not converged
 ST 1-5  9.940e-01  !2147483647  inf  inf  nan!  !1000 0 81!  WARNING: The linear system solution has not converged
 ST 1-6  1.234e+00  [0  nan  nan  nan]  !1000 0 81!  WARNING: The linear system solution has not converged

Reproduction

Run a simulation using the struct equation with the nHK material model fails for 2D triangle elements.

Expected behavior

A simulation using the struct equation with the nHK material model fails for 2D triangle elements should not fail.

Additional context

No response

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[ktbolt]

There are no CI tests for this.

[hanzhao2020]

In sv_struct.cpp function struct_2d line 455:

svMultiPhysics/Code/Source/solver/sv_struct.cpp

Line 455 in fa058d0

lR(0,a) = lR(0,a) + w*(N(a)*ud(0) + Nx(0,a)*P(0,0) + Nx(2,a)*P(0,1));

I think it should be
lR(0,a) = lR(0,a) + w*(N(a)*ud(0) + Nx(0,a)*P(0,0) + Nx(1,a)*P(0,1));

And in mat_models.cpp function compute_pk2cc line 895:

svMultiPhysics/Code/Source/solver/mat_models.cpp

Line 895 in fa058d0

mat_fun::copy_Dm(Dm_2D, Dm, 4, 4);

It looks Dm should have a size of (3,3) for 2D
mat_fun::copy_Dm(Dm_2D, Dm, 3, 3);

Can you try making these two changes and see if the simulation converges?

[ktbolt]

I've found three indexing bugs.

1 read_bc lBc.eDrn is allocated for nsd which is correct. However, in the solver XML file I had <Effective_direction> (1, 0, 0) </Effective_direction> so lBc.eDrn was indexing using effective_direction.size() = 3.

I will add a check that nsd == effective_direction.size()

2 struct_2d An indexing error when computing the local residual lR(0,a).

This was an off-by-one mistake in the original translation from Fortran to C++ made by yours truly.

3 compute_pk2cc The integer arguments passed to mat_fun::copy_Dm(Dm_2D, Dm, 4, 4) are not correct.

[ktbolt]

@hanzhao2020 Thank you for your help on this; found these same problems yesterday !

[ktbolt]

Note that for the mat_fun::copy_Dm(Dm_2D, Dm, 4, 4) call the integer arguments are not only superfluous but are a source of potential mistakes; the Dm_2D and Dm know their sizes so this function can be written as

    void copy_Dm(const MatrixType& mat, Array<double>& dest) 
    {
        if ((mat.rows() != dest.nrows()) || (mat.cols() != dest.ncols())) {
          std::string mat_dims = "(" + std::to_string(mat.rows()) + "x" + std::to_string(mat.cols()) + ")";
          std::string dest_dims = "(" + std::to_string(dest.nrows()) + "x" + std::to_string(dest.ncols()) + ")";
          throw std::runtime_error("[copy_Dm] The 'mat" + mat_dims + "' and 'dest" + dest_dims + "' arrays have incompatible sizes.");
        }

        for (int i = 0; i < mat.rows(); ++i) {
            for (int j = 0; j < mat.cols(); ++j) {
                dest(i, j) = mat(i, j);
            }
        }
    }

[ktbolt]

And in compute_pk2cc the material stiffness tensor Dm is a 3x3 array for 2D but the Eigen Dm_2D array is defined as a 4x4 array

Eigen::Matrix4d Dm_2D = Eigen::Matrix4d::Zero();

[ktbolt]

I've modified compute_pk2cc to set the correct Dm Matrix sizeswrt nsd

const Eigen::Matrix<double, nsd, Eigen::Dynamic> fl, const double ya, Matrix<nsd>& S, Matrix<3*(nsd-1)>& Dm, double& Ja)

2D struct simulations no longer produce NaNs but do not converge

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 ST 1-1  2.939e+00  [0 1.000e+00 1.000e+00 1.706e-04]  !10000 0 96!  WARNING: The linear system solution has not converged
 ST 1-2  5.902e+00  [-41 8.126e-03 8.126e-03 2.023e-02]  !10000 0 96!  WARNING: The linear system solution has not converged
 ST 1-3  8.881e+00  [-13 2.013e-01 2.013e-01 1.821e-03]  !10000 0 96!  WARNING: The linear system solution has not converged
 ST 1-4  1.186e+01  [18 8.445e+00 8.445e+00 2.482e-01]  !10000 0 96!  WARNING: The linear system solution has not converged
 ST 1-5  1.486e+01  !75 5.732e+03 5.732e+03 1.121e-02!  !10000 0 96!  WARNING: The linear system solution has not converged

[ktbolt]

I've found a bug in b_struct_2d computing af; fixed that and now the simulations look good

---------------------------------------------------------------------
 Eq     N-i     T       dB  Ri/R1   Ri/R0    R/Ri     lsIt   dB  %t
---------------------------------------------------------------------
 ST 1-1  2.570e-01  [0 1.000e+00 1.000e+00 9.872e-05]  [508 0 57]
 ST 1-2  6.220e-01  [-67 4.213e-04 4.213e-04 5.805e-05]  [877 -8 70]
 ST 1-3  9.670e-01  [-138 1.130e-07 1.130e-07 9.587e-05]  [815 -1 68]
 ST 1-4  1.142e+00  [-219 1.117e-11 1.117e-11 1.501e-02]  [221 0 37]
 ST 2-1  1.397e+00  [0 1.000e+00 1.301e+00 9.144e-05]  [503 -4 56]
 ST 2-2  1.737e+00  [-65 5.423e-04 7.057e-04 9.855e-05]  [792 0 67]
 ST 2-3  2.031e+00  [-136 1.414e-07 1.840e-07 8.400e-05]  [633 -3 62]
 ST 2-4  2.231e+00  [-218 1.208e-11 1.571e-11 9.684e-03]  [306 -4 44]