SimVascular / SimVascular/svMultiPhysics

Stress output is a zero array in linear elasticity calculation

Open Beginner friendly
#477 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

bug
Dominant language
C++
Stars
45
Forks
60
Avg merge
5d 23h
Merged PRs (30d)
11

Description

Description

The stress output array for linear elasticity (lElas) calculations is consistently a zero array, regardless of non-zero displacement or strain.
The array S (Array<double> S(nsd, nsd)) is not updated within the if (cPhys == EquationType::phys_lElas) block.
The following code is from post.cpp.

        case OutputNameType::outGrp_stress:
        case OutputNameType::outGrp_cauchy: 
        case OutputNameType::outGrp_mises:
          Array<double> sigma(nsd,nsd);
          Array<double> S(nsd,nsd);

          if (cPhys == EquationType::phys_lElas) {
            if (nsd == 3) {
              double detF = lambda*(ed(0) + ed(1) + ed(2));
              sigma(0,0) = detF + 2.0*mu*ed(0);
              sigma(1,1) = detF + 2.0*mu*ed(1);
              sigma(2,2) = detF + 2.0*mu*ed(2);

              sigma(0,1) = mu*ed(3);
              sigma(1,2) = mu*ed(4);
              sigma(2,0) = mu*ed(5);

              sigma(1,0) = sigma(0,1);
              sigma(2,1) = sigma(1,2);
              sigma(0,2) = sigma(2,0);
            } else {  
              double detF = lambda*(ed(0) + ed(1));
              sigma(0,0) = detF + 2.0*mu*ed(0);
              sigma(1,1) = detF + 2.0*mu*ed(1);
              sigma(0,1) = mu*ed(2);
              sigma(1,0) = sigma(0,1);
            }
          } else if (cPhys == EquationType::phys_ustruct) { 
          } else if (cPhys == EquationType::phys_struct) {
          }
          // 2nd Piola-Kirchhoff stress tensor
          if (outGrp == OutputNameType::outGrp_stress) {
            if (nsd == 3) {
              resl(0) = S(0,0);
              resl(1) = S(1,1);
              resl(2) = S(2,2);
              resl(3) = S(0,1);
              resl(4) = S(1,2);
              resl(5) = S(2,0);
            } else { 
              resl(0) = S(0,0);
              resl(1) = S(1,1);
              resl(2) = S(0,1);
            }
Reproduction

The bug can be reproduced by running the test case for linear elasticity.

Expected behavior

If the strain is non-zero, the stress must also be non-zero.

Additional context

No response

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

Contributor guide

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start in post.cpp at the outGrp_stress, outGrp_cauchy, and outGrp_mises handling shown in the issue, then run the linear-elasticity test case. Trace the stress output from the computed strain and verify that non-zero strain produces non-zero stress values in the expected output array.

Written by the indexing model from the issue text.

Assessment

Tech stack
cpp
Domain
backend
Issue type
Bug
Difficulty
2/5
Estimated time
1-3 hours
Activity status
Stale
Clarity
Clearly specified
Newbie friendliness
68/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.