delfem2/test_cpp/fem_stvk_test.cpp at master · nobuyuki83/delfem2 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
//
// Created by Nobuyuki Umetani on 2021/11/20.
//


#include <random>

#include "gtest/gtest.h" // need to be defiend in the beginning

#include "delfem2/vec2.h"
#include "delfem2/fem_stvk.h"

bool RandomTri2(double P[3][2]) {
  namespace dfm2 = delfem2;
  std::mt19937 randomEng(std::random_device{}());
  std::uniform_real_distribution<double> dist_m5p5(-1, 1);
  for (int ip = 0; ip < 3; ++ip) {
    for (int idim = 0; idim < 2; ++idim) {
      P[ip][idim] = dist_m5p5(randomEng);
    }
  }
  if (dfm2::Distance2(P[0], P[1]) < 0.1) { return false; }
  if (dfm2::Distance2(P[1], P[2]) < 0.1) { return false; }
  if (dfm2::Distance2(P[0], P[2]) < 0.1) { return false; }
  double a0 = dfm2::Area_Tri2(P[0], P[1], P[2]);
  if (fabs(a0) < 0.2) return false;
  return true;
}

TEST(fem_stvk, Check_CdC_TriStrain) {
  namespace dfm2 = delfem2;
  std::mt19937 randomEng(std::random_device{}());
  std::uniform_real_distribution<double> dist_m5p5(-5, 5);
// -----
  for (int itr = 0; itr < 200; ++itr) {
    double P[3][2];
    if (!RandomTri2(P)) { continue; }
    const double p[3][3] = {
      {dist_m5p5(randomEng), dist_m5p5(randomEng), dist_m5p5(randomEng)},
      {dist_m5p5(randomEng), dist_m5p5(randomEng), dist_m5p5(randomEng)},
      {dist_m5p5(randomEng), dist_m5p5(randomEng), dist_m5p5(randomEng)}};
    const double eps = 1.0e-5;
    double C[3], dCdp[3][9];
    dfm2::CdC_StVK(C, dCdp, P, p);
    for (int ine = 0; ine < 3; ++ine) {
      for (int idim = 0; idim < 3; ++idim) {
        double p1[3][3];
        for (int i = 0; i < 9; ++i) { (&p1[0][0])[i] = (&p[0][0])[i]; }
        p1[ine][idim] += eps;
        double C1[3], dCdp1[3][9];
        dfm2::CdC_StVK(C1, dCdp1, P, p1);
        EXPECT_NEAR((C1[0] - C[0]) / eps, dCdp[0][ine * 3 + idim], 1.0e-2);
        EXPECT_NEAR((C1[1] - C[1]) / eps, dCdp[1][ine * 3 + idim], 1.0e-2);
        EXPECT_NEAR((C1[2] - C[2]) / eps, dCdp[2][ine * 3 + idim], 1.0e-2);
      }
    }
  }
}

TEST(fem_stvk, Check_CdC_EnergyStVK) {
  namespace dfm2 = delfem2;
  std::mt19937 randomEng(std::random_device{}());
  std::uniform_real_distribution<double> dist_01(0, 1);
  std::uniform_real_distribution<double> dist_12(1, 2);
  //
  for (int itr = 0; itr < 1000; ++itr) {
    double P[3][2];
    if (!RandomTri2(P)) { continue; }
    const double p[3][3] = { //  deformed triangle vertex positions)
      {dist_01(randomEng), dist_01(randomEng), dist_01(randomEng)},
      {dist_01(randomEng), dist_01(randomEng), dist_01(randomEng)},
      {dist_01(randomEng), dist_01(randomEng), dist_01(randomEng)}};
    const double lambda = dist_12(randomEng);
    const double myu = dist_12(randomEng);
    // ---------------------
    double C, dCdp[9];
    dfm2::CdC_EnergyStVK(C, dCdp, P, p, lambda, myu);
    for (int ine = 0; ine < 3; ++ine) {
      for (int idim = 0; idim < 3; ++idim) {
        double eps = 1.0e-8;
        double p1[3][3];
        for (int i = 0; i < 9; ++i) { (&p1[0][0])[i] = (&p[0][0])[i]; }
        p1[ine][idim] += eps;
        double C1, dCdp1[9];
        dfm2::CdC_EnergyStVK(C1, dCdp1, P, p1, lambda, myu);
        EXPECT_NEAR(
          (C1 - C) / eps,
          dCdp[ine * 3 + idim],
          1.0e-3 * (1 + fabs(dCdp[ine * 3 + idim])));
      }
    }
  }
}