-
Notifications
You must be signed in to change notification settings - Fork 1
/
UniformSingleLayer2d_LocalTrigger.cpp
74 lines (61 loc) · 3.26 KB
/
UniformSingleLayer2d_LocalTrigger.cpp
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
#include <catch2/catch.hpp>
#include <xtensor/xrandom.hpp>
#include <FrictionQPotFEM/UniformSingleLayer2d.h>
TEST_CASE("FrictionQPotFEM::UniformSingleLayer2d_LocalTrigger", "UniformSingleLayer2d.h")
{
SECTION("LocalTrigger")
{
size_t N = 6;
double h = xt::numeric_constants<double>::PI;
double L = h * static_cast<double>(N);
GooseFEM::Mesh::Quad4::FineLayer mesh(N, N, h);
auto coor = mesh.coor();
auto conn = mesh.conn();
auto dofs = mesh.dofs();
xt::xtensor<size_t, 1> plastic = mesh.elementsMiddleLayer();
xt::xtensor<size_t, 1> elastic = xt::setdiff1d(xt::arange(mesh.nelem()), plastic);
auto left = mesh.nodesLeftOpenEdge();
auto right = mesh.nodesRightOpenEdge();
xt::view(dofs, xt::keep(right), 0) = xt::view(dofs, xt::keep(left), 0);
xt::view(dofs, xt::keep(right), 1) = xt::view(dofs, xt::keep(left), 1);
auto top = mesh.nodesTopEdge();
auto bottom = mesh.nodesBottomEdge();
size_t nfix = top.size();
xt::xtensor<size_t, 1> iip = xt::empty<size_t>({2 * mesh.ndim() * nfix});
xt::view(iip, xt::range(0 * nfix, 1 * nfix)) = xt::view(dofs, xt::keep(bottom), 0);
xt::view(iip, xt::range(1 * nfix, 2 * nfix)) = xt::view(dofs, xt::keep(bottom), 1);
xt::view(iip, xt::range(2 * nfix, 3 * nfix)) = xt::view(dofs, xt::keep(top), 0);
xt::view(iip, xt::range(3 * nfix, 4 * nfix)) = xt::view(dofs, xt::keep(top), 1);
double c = 1.0;
double G = 1.0;
double K = 10.0 * G;
double rho = G / std::pow(c, 2.0);
double qL = 2.0 * xt::numeric_constants<double>::PI / L;
double qh = 2.0 * xt::numeric_constants<double>::PI / h;
double alpha = std::sqrt(2.0) * qL * c * rho;
double dt = 1.0 / (c * qh) / 10.0;
xt::xtensor<double, 2> epsy = xt::ones<double>(std::array<size_t, 2>{N, 1000});
epsy = xt::cumsum(epsy, 1);
// Initialise system
FrictionQPotFEM::UniformSingleLayer2d::System sys(coor, conn, dofs, iip, elastic, plastic);
sys.setMassMatrix(rho * xt::ones<double>({mesh.nelem()}));
sys.setDampingMatrix(alpha * xt::ones<double>({mesh.nelem()}));
sys.setElastic(K * xt::ones<double>({elastic.size()}), G * xt::ones<double>({elastic.size()}));
sys.setPlastic(K * xt::ones<double>({plastic.size()}), G * xt::ones<double>({plastic.size()}), epsy);
sys.setDt(dt);
{
FrictionQPotFEM::UniformSingleLayer2d::LocalTriggerFineLayer trigger(sys);
xt::xtensor<double, 2> epsy = xt::ones<double>({plastic.size(), size_t(4)});
trigger.setState(sys.Eps(), sys.Sig(), epsy);
xt::xtensor<double, 2> barriers = 5.357607 * xt::ones<double>({plastic.size(), size_t(4)});
REQUIRE(xt::allclose(barriers, trigger.barriers()));
}
{
FrictionQPotFEM::UniformSingleLayer2d::LocalTriggerFineLayerFull trigger(sys);
xt::xtensor<double, 2> epsy = xt::ones<double>({plastic.size(), size_t(4)});
trigger.setState(sys.Eps(), sys.Sig(), epsy);
xt::xtensor<double, 2> barriers = 5.357607 * xt::ones<double>({plastic.size(), size_t(4)});
REQUIRE(xt::allclose(barriers, trigger.barriers()));
}
}
}