PeriDEM 0.3.0
PeriDEM -- Peridynamics-based high-fidelity model for granular media
Loading...
Searching...
No Matches
Namespaces | Functions
main.cpp File Reference
#include "contact/pairForce.h"
#include "contact/policy.h"
#include "contact/wallContact.h"
#include "geom/geomObjects.h"
#include "pd/selfContact.h"
#include "util/io.h"
#include "inp/contactPairDeck.h"
#include <cmath>
#include <cstdlib>
#include <iostream>
Include dependency graph for main.cpp:

Go to the source code of this file.

Namespaces

namespace  anonymous_namespace{main.cpp}
 

Functions

bool anonymous_namespace{main.cpp}::near (double a, double b, double tol=1.e-12)
 
int main ()
 

Function Documentation

◆ main()

int main ( )

Definition at line 31 of file main.cpp.

31 {
32 // Partial volume near contact radius: full at R = Rc - h/2, zero at Rc + h/2.
33 {
34 const double vol = 4.0;
35 const double Rc = 1.0;
36 const double h = 0.2;
37 const double full = contact::correctedContactVolume(vol, 0.5, Rc, h);
38 const double mid = contact::correctedContactVolume(vol, Rc, Rc, h);
39 const double edge = contact::correctedContactVolume(vol, Rc + 0.5 * h, Rc, h);
40 if (!near(full, vol) || !near(mid, 0.5 * vol) || !near(edge, 0.)) {
41 std::cerr << "correctedContactVolume failed: full=" << full
42 << " mid=" << mid << " edge=" << edge << "\n";
43 return 1;
44 }
45 }
46
47 inp::ContactPairDeck deck(/*contactR*/ 1.0, /*computeContactR*/ false,
48 /*dampingOn*/ false, /*frictionOn*/ false,
49 /*Kn*/ 100.0);
50
51 contact::Pair p{deck,
52 util::Point(0., 0., 0.),
53 util::Point(0.5, 0., 0.),
54 util::Point(),
55 util::Point(),
56 0,
57 1,
58 0,
59 1,
60 /*voli*/ 2.0,
61 /*volj*/ 3.0,
62 1.0,
63 1.0,
64 1.e-6,
65 false,
66 false};
67
68 auto law_j = contact::makePairForce("coulomb_simple");
69
70 // Force density spring ∝ Kn*(R-Rc)*Vj (partial volume applied by assembly).
71 const auto fj = law_j->springForce(p);
72 if (!near(fj.d_x, -150.)) {
73 std::cerr << "spring should equal Kn*(R-Rc)*Vj: fj.x=" << fj.d_x << "\n";
74 return 1;
75 }
76
77 // Nodal force would be fj*Vi; density-form integrator uses fj directly.
78 const auto nodal = p.voli * fj;
79 if (!near(nodal.d_x, -300.)) {
80 std::cerr << "nodal force Vi*fj failed: " << nodal.d_x << "\n";
81 return 1;
82 }
83
84 // Node damping fires only when Damping_On; law flags gate COM vs node.
85 {
86 inp::ContactPairDeck dd(/*contactR*/ 1.0, /*computeContactR*/ false,
87 /*dampingOn*/ true, /*frictionOn*/ false,
88 /*Kn*/ 100.0);
89 dd.d_K = 1.e6;
90 dd.d_betan = 0.2;
91 contact::Pair pd{dd,
92 util::Point(0., 0., 0.),
93 util::Point(0.5, 0., 0.),
94 util::Point(1., 0., 0.),
95 util::Point(0., 0., 0.), // approaching
96 0,
97 1,
98 0,
99 1,
100 2.0,
101 3.0,
102 1.0,
103 1.0,
104 1.e-6,
105 false,
106 false};
107 const auto fd = law_j->nodeDampingForce(pd);
108 if (!(fd.length() > 0.)) {
109 std::cerr << "node damping should be nonzero for approaching pair\n";
110 return 1;
111 }
112 if (contact::usesNodeDamping("off") || !contact::usesNodeDamping("node") ||
113 !contact::usesNodeDamping("com_and_node")) {
114 std::cerr << "Damping_Law node flags incorrect\n";
115 return 1;
116 }
117 }
118
119 // Stick-slip: tangential spring capped by mu * |Fn|.
120 {
121 inp::ContactPairDeck fd(/*contactR*/ 1.0, /*computeContactR*/ false,
122 /*dampingOn*/ false, /*frictionOn*/ true,
123 /*Kn*/ 100.0, /*eps*/ 1., /*mu*/ 0.2);
124 contact::Pair q{fd,
125 util::Point(0., 0., 0.),
126 util::Point(0.5, 0., 0.),
127 util::Point(0., 0., 0.),
128 util::Point(0., 10., 0.), // large tangential approach
129 0,
130 1,
131 0,
132 1,
133 1.0,
134 1.0,
135 1.0,
136 1.0,
137 1.e-3,
138 false,
139 false};
140 auto ss = contact::makePairForce("stick_slip");
141 ss->beginStep();
142 // Accumulate enough tangential slip to exceed the Coulomb limit.
143 util::Point fss;
144 for (int k = 0; k < 20; ++k)
145 fss = ss->springForce(q);
146 ss->endStep();
147
148 const double fn = 100.0 * (1.0 - 0.5) * 1.0;
149 const double ft_max = 0.2 * fn;
150 const double ft = std::abs(fss.d_y);
151 if (ft > ft_max + 1.e-9) {
152 std::cerr << "stick-slip exceeded Coulomb cap: ft=" << ft
153 << " ft_max=" << ft_max << "\n";
154 return 1;
155 }
156 if (ft < 0.5 * ft_max) {
157 std::cerr << "stick-slip did not mobilize friction: ft=" << ft << "\n";
158 return 1;
159 }
160 }
161
162 std::cout << "TestContact spring / stick_slip OK\n";
163
164 // Damping_Law selection: COM object vs node-damping flag.
165 {
166 if (contact::makeDamping("off") != nullptr ||
167 contact::makeDamping("node") != nullptr) {
168 std::cerr << "off/node should not create COM damping\n";
169 return 1;
170 }
171 if (contact::makeDamping("com") == nullptr ||
172 contact::makeDamping("com_and_node") == nullptr) {
173 std::cerr << "com/com_and_node should create COM damping\n";
174 return 1;
175 }
176 if (!contact::usesNodeDamping("com_and_node") ||
177 !contact::usesNodeDamping("node") ||
178 contact::usesNodeDamping("com") || contact::usesNodeDamping("off") ||
179 !contact::usesComDamping("com") ||
180 contact::usesComDamping("node")) {
181 std::cerr << "damping law flags incorrect\n";
182 return 1;
183 }
184 }
185
186 std::cout << "TestContact damping laws OK\n";
187
188 // Reference-aware self-contact formula vs broken_bond_kn (natural R = Rc).
189 // Force: Kn * volj * capped_gap / R * yji, gap = R - natural_R, cap = -0.25*natural_R.
190 {
191 auto bb = pd::makeSelfContact("broken_bond_kn");
192 auto rg = pd::makeSelfContact("reference_gap");
193 const util::Point yji(0.4, 0., 0.); // R = 0.4
194 const double volj = 1.0;
195 const double Kn = 10.0;
196 const double Rc = 1.0;
197 const double r0 = 0.5;
198 const auto f_bb = bb->force(yji, volj, Kn, Rc, r0);
199 const auto f_rg = rg->force(yji, volj, Kn, Rc, r0);
200 if (!near(f_bb.d_x, -2.5)) {
201 std::cerr << "broken_bond_kn unexpected: " << f_bb.d_x << "\n";
202 return 1;
203 }
204 if (!near(f_rg.d_x, -1.0)) {
205 std::cerr << "reference_gap unexpected: " << f_rg.d_x << "\n";
206 return 1;
207 }
208 // No force when stretched beyond r0 for reference_gap.
209 const auto f_rg2 = rg->force(util::Point(0.8, 0., 0.), volj, Kn, Rc, r0);
210 if (!near(f_rg2.d_x, 0.) || !near(f_rg2.length(), 0.)) {
211 std::cerr << "reference_gap should be zero when R > r0\n";
212 return 1;
213 }
214 }
215
216 std::cout << "TestContact self-contact OK\n";
217
218 // Wall_Contact factory + geom wallContactQuery (plane / rectangle).
219 {
220 auto meshed = contact::makeWallContact("meshed");
221 auto anal = contact::makeWallContact("analytical_plane");
222 if (!meshed || meshed->skipsMeshedGrainWall() || !anal ||
223 !anal->skipsMeshedGrainWall()) {
224 std::cerr << "Wall_Contact factory flags incorrect\n";
225 return 1;
226 }
227
228 geom::Plane plane(util::Point(0., 1., 0.), util::Point(0., 0., 0.));
229 geom::WallContactHit hit;
230 if (!plane.wallContactQuery(util::Point(0., 0.3, 0.), hit) ||
231 !near(hit.signed_gap, 0.3) || !near(hit.outward_n.d_y, 1.)) {
232 std::cerr << "Plane wallContactQuery free-space failed\n";
233 return 1;
234 }
235 if (!plane.wallContactQuery(util::Point(0., -0.2, 0.), hit) ||
236 !near(hit.signed_gap, -0.2) || !near(hit.outward_n.d_y, 1.)) {
237 std::cerr << "Plane wallContactQuery penetration failed\n";
238 return 1;
239 }
240
241 geom::Rectangle rect(2.0, 1.0, util::Point(0., 0., 0.));
242 if (!rect.wallContactQuery(util::Point(0., 0.8, 0.), hit) ||
243 !near(hit.signed_gap, 0.3) || !near(hit.outward_n.d_y, 1.)) {
244 std::cerr << "Rectangle wallContactQuery failed: gap=" << hit.signed_gap
245 << " ny=" << hit.outward_n.d_y << "\n";
246 return 1;
247 }
248
249 // Analytical wall spring density matches Kn*(gap-Rc)*voli * outward.
250 const double Kn = 100., Rc = 1., vol = 2.;
251 plane.wallContactQuery(util::Point(0., 0.3, 0.), hit);
252 const double gap = hit.signed_gap;
253 auto scalar = Kn * (gap - Rc) * vol;
254 if (scalar > 0.)
255 scalar = 0.;
256 const util::Point f = scalar * (-1. * hit.outward_n);
257 if (!near(f.d_y, 140.) || !near(f.d_x, 0.)) {
258 std::cerr << "analytical wall force failed: f=(" << f.d_x << ","
259 << f.d_y << ")\n";
260 return 1;
261 }
262 }
263
264 std::cout << "TestContact wall contact OK\n";
265
266 // coulomb_simple vs stick_slip differ in tangential response.
267 {
268 inp::ContactPairDeck fd(/*contactR*/ 1.0, /*computeContactR*/ false,
269 /*dampingOn*/ false, /*frictionOn*/ true,
270 /*Kn*/ 100.0, /*eps*/ 1., /*mu*/ 0.2);
271 contact::Pair q{fd,
272 util::Point(0., 0., 0.),
273 util::Point(0.5, 0., 0.),
274 util::Point(0., 0., 0.),
275 util::Point(0., 5., 0.),
276 0,
277 1,
278 0,
279 1,
280 1.0,
281 1.0,
282 1.0,
283 1.0,
284 1.e-3,
285 false,
286 false};
287 auto cs = contact::makePairForce("coulomb_simple");
288 auto ss = contact::makePairForce("stick_slip");
289 const auto f_cs = cs->springForce(q);
290 ss->beginStep();
291 util::Point f_ss;
292 for (int k = 0; k < 5; ++k)
293 f_ss = ss->springForce(q);
294 ss->endStep();
295 // coulomb_simple: instantaneous mu*|Fn| along et; stick_slip builds spring.
296 if (near(f_cs.d_y, f_ss.d_y, 1.e-6) && near(f_cs.d_x, f_ss.d_x, 1.e-6)) {
297 std::cerr << "coulomb_simple and stick_slip should differ in Ft\n";
298 return 1;
299 }
300 }
301
302 std::cout << "TestContact friction laws differ OK\n";
303 return 0;
304}
geom::Plane
Infinite plane: free space is the half-space in the normal direction.
Definition geomObjects.h:428
geom::Rectangle
Defines Rectangle.
Definition geomObjects.h:1254
anonymous_namespace{main.cpp}::near
bool near(double a, double b, double tol=1.e-12)
Definition main.cpp:25
contact::makeDamping
std::unique_ptr< Damping > makeDamping(const std::string &name)
Build COM damping object, or nullptr when the law has no COM term.
Definition policy.cpp:29
contact::makePairForce
std::unique_ptr< PairForce > makePairForce(const std::string &friction_law)
Build pair-force object from friction-law name.
Definition policy.cpp:17
contact::makeWallContact
std::unique_ptr< WallContact > makeWallContact(const std::string &name)
Definition wallContact.cpp:113
contact::correctedContactVolume
double correctedContactVolume(double volj, double Rji, double Rc, double h)
Definition pairForce.cpp:19
contact::usesNodeDamping
bool usesNodeDamping(const std::string &damping_law)
True for node / com_and_node.
Definition policy.cpp:43
contact::usesComDamping
bool usesComDamping(const std::string &damping_law)
True for com / com_and_node.
Definition policy.cpp:39
pd
Definition pdForce.h:18
pd::makeSelfContact
std::unique_ptr< SelfContact > makeSelfContact(const std::string &name)
Build self-contact law named by Model.Self_Contact.
Definition selfContact.cpp:68
contact::Pair
One node-node contact pair. Assembly fills this; the law uses it.
Definition pairForce.h:27
geom::WallContactHit
Result of a wall-contact query against a geom (analytical walls).
Definition geomObjects.h:112
geom::WallContactHit::outward_n
util::Point outward_n
Definition geomObjects.h:115
inp::ContactPairDeck
Structure to read and store particle-particle contact related input data.
Definition contactPairDeck.h:26
util::Point
A structure to represent 3d vectors.
Definition point.h:30
util::Point::d_y
double d_y
the y coordinate
Definition point.h:36
util::Point::d_x
double d_x
the x coordinate
Definition point.h:33

References contact::correctedContactVolume(), inp::ContactPairDeck::d_betan, inp::ContactPairDeck::d_K, util::Point::d_x, util::Point::d_y, contact::makeDamping(), contact::makePairForce(), pd::makeSelfContact(), contact::makeWallContact(), anonymous_namespace{main.cpp}::near(), geom::WallContactHit::outward_n, geom::WallContactHit::signed_gap, contact::usesComDamping(), contact::usesNodeDamping(), geom::Plane::wallContactQuery(), and geom::Rectangle::wallContactQuery().

Here is the call graph for this function: