Skip to content

Commit fa14f76

Browse files
mrp089mrp089
committed
Remove finite differences
1 parent ffcc06a commit fa14f76

File tree

5 files changed

+5
-419
lines changed

5 files changed

+5
-419
lines changed

‎Code/Source/svFSI/ComMod.h

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -935,9 +935,6 @@ class mshType
935935
/// @brief The connectivity array mapping eNoN,nEl --> nNo
936936
Array<int> IEN;
937937

938-
/// @brief map: node -> element
939-
std::vector<std::map<int, std::set<int>>> map_node_ele;
940-
941938
/// @brief gIEN mapper from old to new
942939
Vector<int> otnIEN;
943940

‎Code/Source/svFSI/eq_assem.cpp

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -370,7 +370,6 @@ void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, CmMod& cm_mod, const mshT
370370

371371
case EquationType::phys_gr:
372372
gr::construct_gr(com_mod, lM, Dg, false);
373-
// gr::construct_gr_fd_global(com_mod, lM, Dg);
374373
break;
375374

376375
case EquationType::phys_ustruct:

‎Code/Source/svFSI/gr_struct.cpp

Lines changed: 5 additions & 242 deletions
Original file line numberDiff line numberDiff line change
@@ -78,254 +78,17 @@ void construct_gr(ComMod &com_mod, const mshType &lM, const Array<double> &Dg,
7878
lR = 0.0;
7979
lK = 0.0;
8080

81-
// Evaluate solid equations
82-
if (eval_fd) {
83-
eval_gr_fd_ele(e, com_mod, lM, e_Dg, ptr, lR, lK);
84-
} else {
85-
// Update G&R internal variables
86-
eval_gr(e, com_mod, lM, Dg, ptr, lR, lK, false, false);
87-
88-
// Compute stress and tangent
89-
eval_gr(e, com_mod, lM, Dg, ptr, lR, lK, true, true);
90-
}
81+
// Update G&R internal variables
82+
eval_gr(e, com_mod, lM, Dg, ptr, lR, lK, false, false);
83+
84+
// Compute stress and tangent
85+
eval_gr(e, com_mod, lM, Dg, ptr, lR, lK, true, true);
9186

9287
// Assemble into global residual and tangent
9388
lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR);
9489
}
9590
}
9691

97-
/// @brief Finite difference on each element
98-
void eval_gr_fd_ele(const int &e, ComMod &com_mod, const mshType &lM,
99-
Array<double> &Dg, Vector<int> &ptr, Array<double> &lR,
100-
Array3<double> &lK) {
101-
// Get dimensions
102-
const int eNoN = lM.eNoN;
103-
const int dof = com_mod.dof;
104-
const int tDof = com_mod.tDof;
105-
const int tnNo = com_mod.tnNo;
106-
107-
// Set step size for finite difference
108-
const double eps = 1.0e-10;
109-
110-
// Time integration parameters
111-
int Ac;
112-
int cEq = com_mod.cEq;
113-
auto &eq = com_mod.eq[cEq];
114-
115-
// Initialize residual and tangent
116-
Array<double> lRp(dof, eNoN);
117-
Array<double> dlR(dof, eNoN);
118-
Array3<double> lK_dummy(dof * dof, eNoN, eNoN);
119-
lRp = 0.0;
120-
lK_dummy = 0.0;
121-
122-
// Central evaluation
123-
eval_gr(e, com_mod, lM, Dg, ptr, lR, lK_dummy, true, false);
124-
125-
// Finite differences
126-
for (int i = 0; i < dof; ++i) {
127-
for (int a = 0; a < eNoN; ++a) {
128-
dlR = 0.0;
129-
130-
// Global node ID
131-
Ac = lM.IEN(a, e);
132-
133-
// Perturb
134-
Dg(i, Ac) += eps;
135-
136-
// Displacement
137-
lRp = 0.0;
138-
eval_gr(e, com_mod, lM, Dg, ptr, lRp, lK_dummy, true, false);
139-
dlR += (lRp - lR) / eps;
140-
141-
// Restore
142-
Dg(i, Ac) -= eps;
143-
144-
// Assign to tangent matrix
145-
for (int j = 0; j < dof; ++j) {
146-
for (int b = 0; b < eNoN; ++b) {
147-
lK(j * dof + i, b, a) = dlR(j, b);
148-
}
149-
}
150-
}
151-
}
152-
}
153-
154-
void construct_gr_fd_global(ComMod &com_mod, const mshType &lM,
155-
const Array<double> &Dg) {
156-
// Get dimensions
157-
const int eNoN = lM.eNoN;
158-
const int dof = com_mod.dof;
159-
const int tDof = com_mod.tDof;
160-
const int tnNo = com_mod.tnNo;
161-
162-
// Set step size for finite difference
163-
const double eps = 1.0e-8;
164-
165-
// Make editable copy
166-
Array<double> e_Dg(tDof, tnNo);
167-
e_Dg = Dg;
168-
169-
// Initialize residual and tangent
170-
Vector<int> ptr(eNoN);
171-
Array<double> lR(dof, eNoN);
172-
Array3<double> lK(dof * dof, eNoN, eNoN);
173-
174-
// Residual evaluation
175-
eval_gr_fd_global(com_mod, lM, Dg, 1.0 / eps);
176-
177-
// Loop global nodes
178-
for (int Ac = 0; Ac < tnNo; ++Ac) {
179-
// Central evaluation
180-
eval_gr_fd_global(com_mod, lM, Dg, 1.0 / eps, Ac);
181-
182-
// Loop DOFs
183-
for (int i = 0; i < dof; ++i) {
184-
// Perturb solution vectors
185-
e_Dg(i, Ac) += eps;
186-
187-
// Perturbed evaluations
188-
eval_gr_fd_global(com_mod, lM, e_Dg, 1.0 / eps, Ac, i);
189-
190-
// Restore solution vectors
191-
e_Dg(i, Ac) = Dg(i, Ac);
192-
}
193-
}
194-
}
195-
196-
void eval_gr_fd_global(ComMod &com_mod, const mshType &lM,
197-
const Array<double> &Dg, const double eps, const int dAc,
198-
const int dj) {
199-
using namespace consts;
200-
201-
// Check if the residual should be assembled
202-
const bool residual = (dAc == -1) && (dj == -1);
203-
204-
// Check if this is the central evaluation
205-
const bool central = (dAc != -1) && (dj == -1);
206-
207-
// Get dimensions
208-
const int eNoN = lM.eNoN;
209-
const int dof = com_mod.dof;
210-
211-
// Initialize residual and tangent
212-
Vector<int> ptr_dummy(eNoN);
213-
Array<double> lR_dummy(dof, eNoN);
214-
Array3<double> lK_dummy(dof * dof, eNoN, eNoN);
215-
ptr_dummy = 0;
216-
lR_dummy = 0.0;
217-
lK_dummy = 0.0;
218-
219-
// Select element set to evaluate
220-
std::set<int> elements;
221-
if (residual) {
222-
// Residual evaluation: evaluate all elements
223-
for (int i = 0; i < lM.nEl; ++i) {
224-
elements.insert(i);
225-
}
226-
} else {
227-
// Pick elements according to smoothing algorithm
228-
elements = lM.map_node_ele[1].at(dAc);
229-
}
230-
231-
// Update internal G&R variables without evaluating stress or tangent
232-
for (int e : elements) {
233-
eval_gr(e, com_mod, lM, Dg, ptr_dummy, lR_dummy, lK_dummy, false, false);
234-
}
235-
236-
// Index of collagen mass fraction multiplier
237-
const int igr = 37;
238-
239-
// Initialize arrays
240-
Vector<double> grInt_a(lM.gnNo);
241-
Vector<double> grInt_n(lM.gnNo);
242-
grInt_a = 0.0;
243-
grInt_n = 0.0;
244-
245-
int Ac;
246-
double w;
247-
double val;
248-
Vector<double> N(lM.eNoN);
249-
250-
// Project: integration point -> nodes
251-
for (int g = 0; g < lM.nG; g++) {
252-
w = lM.w(g);
253-
N = lM.N.col(g);
254-
for (int e : elements) {
255-
val = com_mod.grInt(e, g, igr);
256-
for (int a = 0; a < lM.eNoN; a++) {
257-
Ac = lM.IEN(a, e);
258-
// todo: add jacobian
259-
grInt_n(Ac) += w * N(a) * val;
260-
grInt_a(Ac) += w * N(a);
261-
}
262-
}
263-
}
264-
265-
// Project: nodes -> integration points
266-
for (int e : elements) {
267-
for (int g = 0; g < lM.nG; g++) {
268-
N = lM.N.col(g);
269-
val = 0.0;
270-
for (int a = 0; a < lM.eNoN; a++) {
271-
Ac = lM.IEN(a, e);
272-
val += N(a) * grInt_n(Ac) / grInt_a(Ac);
273-
}
274-
com_mod.grInt(e, g, igr) = val;
275-
}
276-
}
277-
278-
// Store internal G&R variables
279-
if (residual) {
280-
com_mod.grInt_orig = com_mod.grInt;
281-
}
282-
283-
// Initialzie arrays for Finite Difference (FD)
284-
Vector<int> ptr_row(eNoN);
285-
Vector<int> ptr_col(1);
286-
Array<double> lR(dof, eNoN);
287-
Array3<double> lK(dof * dof, eNoN, 1);
288-
289-
// Assemble only the FD node
290-
ptr_col = dAc;
291-
292-
// Loop over all elements of mesh
293-
for (int e : elements) {
294-
// Reset
295-
ptr_row = 0;
296-
lR = 0.0;
297-
lK = 0.0;
298-
299-
// Evaluate solid equations (with smoothed internal G&R variables)
300-
eval_gr(e, com_mod, lM, Dg, ptr_row, lR, lK_dummy, true, false);
301-
302-
// Assemble into global residual
303-
if (residual) {
304-
lhsa_ns::do_assem_residual(com_mod, lM.eNoN, ptr_row, lR);
305-
continue;
306-
}
307-
308-
// Components of FD: central and difference
309-
for (int a = 0; a < eNoN; ++a) {
310-
for (int i = 0; i < dof; ++i) {
311-
if (central) {
312-
for (int j = 0; j < dof; ++j) {
313-
lK(i * dof + j, a, 0) = -lR(i, a) * eps;
314-
}
315-
} else {
316-
lK(i * dof + dj, a, 0) = lR(i, a) * eps;
317-
}
318-
}
319-
}
320-
321-
// Assemble into global tangent
322-
lhsa_ns::do_assem_tangent(com_mod, lM.eNoN, 1, ptr_row, ptr_col, lK);
323-
}
324-
325-
// Restore internal G&R variables
326-
com_mod.grInt = com_mod.grInt_orig;
327-
}
328-
32992
/// @brief
33093
void eval_gr(const int &e, ComMod &com_mod, const mshType &lM,
33194
const Array<double> &Dg, Vector<int> &ptr, Array<double> &lR,

‎Code/Source/svFSI/gr_struct.h

Lines changed: 0 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -45,21 +45,6 @@ namespace gr {
4545
void construct_gr(ComMod &com_mod, const mshType &lM, const Array<double> &Dg,
4646
bool eval_fd = true);
4747

48-
void construct_gr_fd_global(ComMod &com_mod, const mshType &lM,
49-
const Array<double> &Dg);
50-
51-
void eval_gr_fd_global(ComMod &com_mod, const mshType &lM,
52-
const Array<double> &Dg, const double eps = 0.0,
53-
const int dAc = -1, const int di = -1);
54-
55-
void eval_gr_fd_global_phic(ComMod &com_mod, const mshType &lM,
56-
const Array<double> &Dg, const double eps = 0.0,
57-
const int dAc = -1, const int di = -1);
58-
59-
void eval_gr_fd_ele(const int &e, ComMod &com_mod, const mshType &lM,
60-
Array<double> &Dg, Vector<int> &ptr, Array<double> &lR,
61-
Array3<double> &lK);
62-
6348
void eval_gr(const int &e, ComMod &com_mod, const mshType &lM,
6449
const Array<double> &Dg, Vector<int> &ptr, Array<double> &lR,
6550
Array3<double> &lK, const bool eval_s = true,

0 commit comments

Comments
 (0)