11363842044

Committed 16 Oct 2024 10:30AM UTC coverage: 0.794% (-8.0%) from 8.812%

Build # 11363842044

Build Type

push

github

Committed by

Pablo1990

Commit Message

Big refactor on interface type

might affect performance, but minimise errors

Run Details

0 of 1002 branches covered (0.0%)

Branch coverage included in aggregate %.

1 of 64 new or added lines in 10 files covered. (1.56%)

524 existing lines in 33 files now uncovered.

51 of 5423 relevant lines covered (0.94%)

0.01 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/pyVertexModel/algorithm/newtonRaphson.py

UNCOV 1	import logging	×
2
UNCOV 3	import numpy as np	×
UNCOV 4	from scipy.optimize import minimize, Bounds	×
5
UNCOV 6	from src.pyVertexModel.Kg.kgContractility import KgContractility	×
UNCOV 7	from src.pyVertexModel.Kg.kgContractility_external import KgContractilityExternal	×
UNCOV 8	from src.pyVertexModel.Kg.kgSubstrate import KgSubstrate	×
UNCOV 9	from src.pyVertexModel.Kg.kgSurfaceCellBasedAdhesion import KgSurfaceCellBasedAdhesion	×
UNCOV 10	from src.pyVertexModel.Kg.kgTriAREnergyBarrier import KgTriAREnergyBarrier	×
UNCOV 11	from src.pyVertexModel.Kg.kgTriEnergyBarrier import KgTriEnergyBarrier	×
UNCOV 12	from src.pyVertexModel.Kg.kgViscosity import KgViscosity	×
UNCOV 13	from src.pyVertexModel.Kg.kgVolume import KgVolume	×
14
UNCOV 15	logger = logging.getLogger("pyVertexModel")	×
16
17
UNCOV 18	def solve_remodeling_step(geo_0, geo_n, geo, dofs, c_set):	×
19	"""
20	This function solves local problem to obtain the position of the newly remodeled vertices with prescribed settings
21	(Set.***_LP), e.g. Set.lambda_LP.
22	:param geo_0:
23	:param geo_n:
24	:param geo:
25	:param dofs:
26	:param c_set:
27	:return:
28	"""
29
30	logger.info('=====>> Solving Local Problem....')	×
31	geo.remodelling = True	×
32	original_nu = c_set.nu	×
33	original_nu0 = c_set.nu0	×
34	original_lambdaB = c_set.lambdaB	×
35	original_lambdaR = c_set.lambdaR	×
36
37	c_set.MaxIter = c_set.MaxIter0 * 3	×
38	c_set.lambdaB = original_lambdaB * 2	×
39	c_set.lambdaR = original_lambdaR * 0.0001	×
40
41	for n_id, nu_factor in enumerate(np.linspace(10, 1, 3)):	×
42	c_set.nu0 = original_nu0 * nu_factor	×
43	c_set.nu = original_nu * nu_factor	×
44	g, k, _, _ = KgGlobal(geo_0, geo_n, geo, c_set)	×
45
46	dy = np.zeros(((geo.numY + geo.numF + geo.nCells) * 3, 1), dtype=np.float64)	×
47	dyr = np.linalg.norm(dy[dofs.remodel, 0])	×
48	gr = np.linalg.norm(g[dofs.remodel])	×
49	logger.info(	×
50	f'Local Problem ->Iter: 0, \|\|gr\|\|= {gr:.3e} \|\|dyr\|\|= {dyr:.3e} nu/nu0={c_set.nu / c_set.nu0:.3e} '
51	f'dt/dt0={c_set.dt / c_set.dt0:.3g}')
52
53	geo, g, k, energy, c_set, gr, dyr, dy = newton_raphson(geo_0, geo_n, geo, dofs, c_set, k, g, -1, -1)	×
54
55	if gr > c_set.tol or dyr > c_set.tol or np.any(np.isnan(g[dofs.Free])) or np.any(np.isnan(dy[dofs.Free])):	×
56	logger.info(f'Local Problem did not converge after {c_set.iter} iterations.')	×
57	has_converged = False	×
58	break	×
59	else:
60	logger.info(f'=====>> Local Problem converged in {c_set.iter} iterations.')	×
61	has_converged = True	×
62
63	geo.remodelling = False	×
64
65	c_set.MaxIter = c_set.MaxIter0	×
66	c_set.nu = original_nu	×
67	c_set.nu0 = original_nu0	×
68	c_set.lambdaB = original_lambdaB	×
69	c_set.lambdaR = original_lambdaR	×
70
71	return geo, c_set, has_converged	×
72
73
UNCOV 74	def newton_raphson(Geo_0, Geo_n, Geo, Dofs, Set, K, g, numStep, t, implicit_method=True):	×
75	"""
76	Newton-Raphson method
77	:param Geo_0:
78	:param Geo_n:
79	:param Geo:
80	:param Dofs:
81	:param Set:
82	:param K:
83	:param g:
84	:param numStep:
85	:param t:
86	:return:
87	"""
88	if Geo.remodelling:	×
89	dof = Dofs.remodel	×
90	else:
91	dof = Dofs.Free	×
92
93	dy = np.zeros(((Geo.numY + Geo.numF + Geo.nCells) * 3, 1), dtype=np.float64)	×
94	dyr = np.linalg.norm(dy[dof, 0])	×
95	gr = np.linalg.norm(g[dof])	×
96	gr0 = gr	×
97
98	logger.info(f"Step: {numStep}, Iter: 0 \|\|gr\|\|= {gr} \|\|dyr\|\|= {dyr} dt/dt0={Set.dt / Set.dt0:.3g}")	×
99
100	Set.iter = 1	×
101	auxgr = np.zeros(3, dtype=np.float64)	×
102	auxgr[0] = gr	×
103	ig = 0	×
104	Energy = 0	×
105
106	if implicit_method is True:	×
107	while (gr > Set.tol or dyr > Set.tol) and Set.iter < Set.MaxIter:	×
108	Energy, K, dyr, g, gr, ig, auxgr, dy = newton_raphson_iteration(Dofs, Geo, Geo_0, Geo_n, K, Set, auxgr, dof,	×
109	dy, g, gr0, ig, numStep, t)
110	else:
111	Geo, dy, gr = newton_raphson_iteration_explicit(Geo, Set, dof, dy, g)	×
112	dyr = 0	×
113
114	return Geo, g, K, Energy, Set, gr, dyr, dy	×
115
116
UNCOV 117	def newton_raphson_iteration(Dofs, Geo, Geo_0, Geo_n, K, Set, aux_gr, dof, dy, g, gr0, ig, numStep, t):	×
118	"""
119	Perform a single iteration of the Newton-Raphson method for solving a system of nonlinear equations.
120
121	Parameters:
122	Dofs (object): Object containing information about the degrees of freedom in the system.
123	Geo (object): Object containing the current geometry of the system.
124	Geo_0 (object): Object containing the initial geometry of the system.
125	Geo_n (object): Object containing the geometry of the system at the previous time step.
126	K (ndarray): The Jacobian matrix of the system.
127	Set (object): Object containing various settings for the simulation.
128	aux_gr (ndarray): Array containing the norms of the gradient at the previous three steps.
129	dof (ndarray): Array containing the indices of the free degrees of freedom in the system.
130	dy (ndarray): Array containing the current guess for the solution of the system.
131	g (ndarray): The gradient of the system.
132	gr0 (float): The norm of the gradient at the start of the current time step.
133	ig (int): Index used to cycle through the elements of aux_gr.
134	numStep (int): The current time step number.
135	t (float): The current time.
136
137	Returns:
138	energy_total (float): The total energy of the system after the current iteration.
139	K (ndarray): The updated Jacobian matrix of the system.
140	dyr (float): The norm of the change in the solution guess during the current iteration.
141	g (ndarray): The updated gradient of the system.
142	gr (float): The norm of the updated gradient of the system.
143	ig (int): The updated index used to cycle through the elements of aux_gr.
144	aux_gr (ndarray): The updated array containing the norms of the gradient at the previous three steps.
145	dy (ndarray): The updated guess for the solution of the system.
146	"""
147	dy[dof, 0] = ml_divide(K, dof, g)	×
148
149	alpha = line_search(Geo_0, Geo_n, Geo, dof, Set, g, dy)	×
150	dy_reshaped = np.reshape(dy * alpha, (Geo.numF + Geo.numY + Geo.nCells, 3))	×
151	Geo.update_vertices(dy_reshaped)	×
152	Geo.update_measures()	×
153	g, K, energy_total, _ = KgGlobal(Geo_0, Geo_n, Geo, Set)	×
154
155	dyr = np.linalg.norm(dy[dof, 0])	×
156	gr = np.linalg.norm(g[dof])	×
157	logger.info(f"Step: {numStep}, Iter: {Set.iter}, Time: {t} \|\|gr\|\|= {gr:.3e} \|\|dyr\|\|= {dyr:.3e} alpha= {alpha:.3e}"	×
158	f" nu/nu0={Set.nu / Set.nu0:.3g}")
159	Set.iter += 1	×
160
161	# Checking if the three previous steps are very similar. Thus, the solution is not converging
162	aux_gr[ig] = gr	×
163	ig = 0 if ig == 2 else ig + 1	×
164	if (	×
165	all(aux_gr[i] != 0 for i in range(3)) and
166	abs(aux_gr[0] - aux_gr[1]) / aux_gr[0] < 1e-3 and
167	abs(aux_gr[0] - aux_gr[2]) / aux_gr[0] < 1e-3 and
168	abs(aux_gr[2] - aux_gr[1]) / aux_gr[2] < 1e-3
169	) or (
170	gr0 != 0 and abs((gr0 - gr) / gr0) > 1e3
171	):
172	Set.iter = Set.MaxIter	×
173
174	return energy_total, K, dyr, g, gr, ig, aux_gr, dy	×
175
176
UNCOV 177	def newton_raphson_iteration_explicit(Geo, Set, dof, dy, g):	×
178	"""
179	Explicit update method
180	:param Geo_0:
181	:param Geo_n:
182	:param Geo:
183	:param Dofs:
184	:param Set:
185	:return:
186	"""
187	dy[dof, 0] = -Set.dt / Set.nu * g[dof]	×
188
189	# Update border ghost nodes with the same displacement as the corresponding real nodes
190	dy = map_vertices_periodic_boundaries(Geo, dy)	×
191
192	dy_reshaped = np.reshape(dy, (Geo.numF + Geo.numY + Geo.nCells, 3))	×
193	Geo.update_vertices(dy_reshaped)	×
194	Geo.update_measures()	×
195
196	g, energies = gGlobal(Geo, Geo, Geo, Set, Set.implicit_method)	×
197	gr = np.linalg.norm(g[dof])	×
198	Set.iter = Set.MaxIter	×
199
200	return Geo, dy, gr	×
201
202
UNCOV 203	def map_vertices_periodic_boundaries(Geo, dy):	×
204	"""
205	Update the border ghost nodes with the same displacement as the corresponding real nodes
206	:param Geo:
207	:param dy:
208	:return:
209	"""
210	for numCell in Geo.BorderCells:	×
211	cell = Geo.Cells[numCell]	×
212
213	tets_to_check = cell.T[np.any(np.isin(cell.T, Geo.BorderGhostNodes), axis=1)]	×
214	global_ids_to_update = cell.globalIds[np.any(np.isin(cell.T, Geo.BorderGhostNodes), axis=1)]	×
215
216	for i, tet in enumerate(tets_to_check):	×
217	for j, node in enumerate(tet):	×
218	if node in Geo.BorderGhostNodes:	×
219	global_id = global_ids_to_update[i]	×
220
221	# Check if the object has the attribute opposite_cell
222	if not hasattr(Geo.Cells[int(node)], 'opposite_cell'):	×
223	return dy	×
224
225	# Iterate over each element in tet and access the opposite_cell attribute
226	opposite_cell = Geo.Cells[Geo.Cells[int(node)].opposite_cell]	×
227	opposite_cells = np.unique(	×
228	[Geo.Cells[int(t)].opposite_cell for t in tet if Geo.Cells[int(t)].opposite_cell is not None])
229
230	# Get the most similar tetrahedron in the opposite cell
231	if len(opposite_cells) < 3:	×
232	possible_tets = opposite_cell.T[	×
233	[np.sum(np.isin(opposite_cells, tet)) == len(opposite_cells) for tet in opposite_cell.T]]
234	else:
235	possible_tets = opposite_cell.T[	×
236	[np.sum(np.isin(opposite_cells, tet)) >= 2 for tet in opposite_cell.T]]
237
238	if possible_tets.shape[0] == 0:	×
239	possible_tets = opposite_cell.T[	×
240	[np.sum(np.isin(opposite_cells, tet)) == 1 for tet in opposite_cell.T]]
241
242	if possible_tets.shape[0] > 0:	×
243	# Filter the possible tets to get its correct location (XgTop, XgBottom or XgLateral)
244	if possible_tets.shape[0] > 1:	×
245	old_possible_tets = possible_tets	×
246	scutoid = False	×
247	if np.isin(tet, Geo.XgTop).any():	×
248	possible_tets = possible_tets[np.isin(possible_tets, Geo.XgTop).any(axis=1)]	×
249	elif np.isin(tet, Geo.XgBottom).any():	×
250	possible_tets = possible_tets[np.isin(possible_tets, Geo.XgBottom).any(axis=1)]	×
251	else:
252	# Get the tets that are not in the top or bottom
253	possible_tets = possible_tets[	×
254	np.logical_not(np.isin(possible_tets, Geo.XgTop).any(axis=1))]
255	scutoid = True	×
256
257	if possible_tets.shape[0] == 0:	×
258	possible_tets = old_possible_tets	×
259
260	# Compare tets with the number of ghost nodes
261	if possible_tets.shape[0] > 1 and not scutoid:	×
262	old_possible_tets = possible_tets	×
263	# Get the tets that have the same number of ghost nodes as the original tet
264	possible_tets = possible_tets[	×
265	np.sum(np.isin(tet[tet != node], Geo.XgID)) == np.sum(np.isin(possible_tets, Geo.XgID),
266	axis=1)]
267
268	if possible_tets.shape[0] == 0:	×
269	possible_tets = old_possible_tets	×
270
271	if possible_tets.shape[0] > 1 and not scutoid:	×
272	old_possible_tets = possible_tets	×
273	# Get the tet that has the same number of ghost nodes as the original tet
274	possible_tets = possible_tets[np.sum(np.isin(tet[tet != node], Geo.XgID)) == np.sum(	×
275	np.isin(possible_tets, np.concatenate([Geo.XgTop, Geo.XgBottom])), axis=1)]
276
277	if possible_tets.shape[0] == 0:	×
278	possible_tets = old_possible_tets	×
279	else:
280	print('Error in Tet ID: ', tet)	×
281
282	if possible_tets.shape[0] == 1:	×
283	opposite_global_id = opposite_cell.globalIds[	×
284	np.where(np.all(opposite_cell.T == possible_tets[0], axis=1))[0][0]]
285	dy[global_id * 3:global_id * 3 + 3, 0] = dy[opposite_global_id * 3:opposite_global_id * 3 + 3,	×
286	0]
287	elif possible_tets.shape[0] > 1:	×
288	# TODO: what happens if it is an scutoid
289	if scutoid:	×
290	avg_dy = np.zeros(3)	×
291	# Average of the dys
292	for c_tet in possible_tets:	×
293	opposite_global_id = opposite_cell.globalIds[	×
294	np.where(np.all(opposite_cell.T == c_tet, axis=1))[0][0]]
295	avg_dy += dy[opposite_global_id * 3:opposite_global_id * 3 + 3, 0]	×
296	avg_dy /= possible_tets.shape[0]	×
297
298	# NOW IT IS 0, WITH THE AVERAGE IT THERE WERE ERRORS
299	dy[global_id * 3:global_id * 3 + 3, 0] = 0	×
300	else:
301	opposite_global_id = opposite_cell.globalIds[	×
302	np.where(np.all(opposite_cell.T == possible_tets[0], axis=1))[0][0]]
303	dy[global_id * 3:global_id * 3 + 3, 0] = dy[	×
304	opposite_global_id * 3:opposite_global_id * 3 + 3,
305	0]
306	return dy	×
307
308
UNCOV 309	def ml_divide(K, dof, g):	×
310	"""
311	Solve the linear system K * dy = g
312	:param K:
313	:param dof:
314	:param g:
315	:return:
316	"""
317	# dy[dof] = kg_functions.mldivide_np(K[np.ix_(dof, dof)], g[dof])
318	return -np.linalg.solve(K[np.ix_(dof, dof)], g[dof])	×
319
320
UNCOV 321	def line_search(Geo_0, Geo_n, geo, dof, Set, gc, dy):	×
322	"""
323	Line search to find the best alpha to minimize the energy
324	:param Geo_0: Initial geometry.
325	:param Geo_n: Geometry at the previous time step
326	:param geo: Current geometry
327	:param Dofs: Dofs object
328	:param Set: Set object
329	:param gc: Gradient at the current step
330	:param dy: Displacement at the current step
331	:return: alpha
332	"""
333	dy_reshaped = np.reshape(dy, (geo.numF + geo.numY + geo.nCells, 3))	×
334
335	# Create a copy of geo to not change the original one
336	Geo_copy = geo.copy()	×
337
338	Geo_copy.update_vertices(dy_reshaped)	×
339	Geo_copy.update_measures()	×
340
341	g, _ = gGlobal(Geo_0, Geo_n, Geo_copy, Set)	×
342
343	gr0 = np.linalg.norm(gc[dof])	×
344	gr = np.linalg.norm(g[dof])	×
345
346	if gr0 < gr:	×
347	R0 = np.dot(dy[dof, 0], gc[dof])	×
348	R1 = np.dot(dy[dof, 0], g[dof])	×
349
350	R = R0 / R1	×
351	alpha1 = (R / 2) + np.sqrt((R / 2) ** 2 - R)	×
352	alpha2 = (R / 2) - np.sqrt((R / 2) ** 2 - R)	×
353
354	if np.isreal(alpha1) and 2 > alpha1 > 1e-3:	×
355	alpha = alpha1	×
356	elif np.isreal(alpha2) and 2 > alpha2 > 1e-3:	×
357	alpha = alpha2	×
358	else:
359	alpha = 0.1	×
360	else:
361	alpha = 1	×
362
363	return alpha	×
364
365
UNCOV 366	def KgGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True):	×
367	"""
368	Compute the global Jacobian matrix and the global gradient
369	:param Geo_0:
370	:param Geo_n:
371	:param Geo:
372	:param Set:
373	:param implicit_method:
374	:return:
375	"""
376
377	# Surface Energy
378	kg_SA = KgSurfaceCellBasedAdhesion(Geo)	×
379	kg_SA.compute_work(Geo, Set)	×
380	g = kg_SA.g	×
381	K = kg_SA.K	×
382	energy_total = kg_SA.energy	×
383	energies = {"Surface": kg_SA.energy}	×
384
385	# Volume Energy
386	kg_Vol = KgVolume(Geo)	×
387	kg_Vol.compute_work(Geo, Set)	×
388	K += kg_Vol.K	×
389	g += kg_Vol.g	×
390	energy_total += kg_Vol.energy	×
391	energies["Volume"] = kg_Vol.energy	×
392
393	# # TODO: Plane Elasticity
394	# if Set.InPlaneElasticity:
395	# gt, Kt, EBulk = KgBulk(geo_0, Geo, Set)
396	# K += Kt
397	# g += gt
398	# E += EBulk
399	# Energies["Bulk"] = EBulk
400
401	# Bending Energy
402	# TODO
403
404	# Propulsion Forces
405	# TODO
406
407	# Contractility
408	if Set.Contractility:	×
409	kg_lt = KgContractility(Geo)	×
410	kg_lt.compute_work(Geo, Set)	×
411	g += kg_lt.g	×
412	K += kg_lt.K	×
413	energy_total += kg_lt.energy	×
414	energies["Contractility"] = kg_lt.energy	×
415
416	if Set.Contractility_external:	×
417	kg_ext = KgContractilityExternal(Geo)	×
418	kg_ext.compute_work(Geo, Set)	×
419	g += kg_ext.g	×
420	K += kg_ext.K	×
421	energy_total += kg_ext.energy	×
422	energies["Contractility_external"] = kg_ext.energy	×
423
424	# Substrate
425	if Set.Substrate == 2:	×
426	kg_subs = KgSubstrate(Geo)	×
427	kg_subs.compute_work(Geo, Set)	×
428	g += kg_subs.g	×
429	K += kg_subs.K	×
430	energy_total += kg_subs.energy	×
431	energies["Substrate"] = kg_subs.energy	×
432
433	# Triangle Energy Barrier
434	if Set.EnergyBarrierA:	×
435	kg_Tri = KgTriEnergyBarrier(Geo)	×
436	kg_Tri.compute_work(Geo, Set)	×
437	g += kg_Tri.g	×
438	K += kg_Tri.K	×
439	energy_total += kg_Tri.energy	×
440	energies["TriEnergyBarrier"] = kg_Tri.energy	×
441
442	# Triangle Energy Barrier Aspect Ratio
443	if Set.EnergyBarrierAR:	×
444	kg_TriAR = KgTriAREnergyBarrier(Geo)	×
445	kg_TriAR.compute_work(Geo, Set)	×
446	g += kg_TriAR.g	×
447	K += kg_TriAR.K	×
448	energy_total += kg_TriAR.energy	×
449	energies["TriEnergyBarrierAR"] = kg_TriAR.energy	×
450
451	if implicit_method is True:	×
452	# Viscous Energy
453	kg_Viscosity = KgViscosity(Geo)	×
454	kg_Viscosity.compute_work(Geo, Set, Geo_n)	×
455	g += kg_Viscosity.g	×
456	K += kg_Viscosity.K	×
457	energy_total = kg_Viscosity.energy	×
458	energies["Viscosity"] = kg_Viscosity.energy	×
459
460	return g, K, energy_total, energies	×
461
462
UNCOV 463	def gGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True):	×
464	"""
465	Compute the global gradient
466	:param Geo_0:
467	:param Geo_n:
468	:param Geo:
469	:param Set:
470	:param implicit_method:
471	:return:
472	"""
473	g = np.zeros((Geo.numY + Geo.numF + Geo.nCells) * 3, dtype=np.float64)	×
474	energies = {}	×
475
476	# Surface Energy
477	if Set.lambdaS1 > 0:	×
478	kg_SA = KgSurfaceCellBasedAdhesion(Geo)	×
479	kg_SA.compute_work(Geo, Set, None, False)	×
480	g += kg_SA.g	×
481	energies["Surface"] = kg_SA.energy	×
482
483	# Volume Energy
484	if Set.lambdaV > 0:	×
485	kg_Vol = KgVolume(Geo)	×
486	kg_Vol.compute_work(Geo, Set, None, False)	×
487	g += kg_Vol.g[:]	×
488	energies["Volume"] = kg_Vol.energy	×
489
490	if implicit_method is True:	×
491	# Viscous Energy
492	kg_Viscosity = KgViscosity(Geo)	×
493	kg_Viscosity.compute_work(Geo, Set, Geo_n, False)	×
494	g += kg_Viscosity.g	×
495	energies["Viscosity"] = kg_Viscosity.energy	×
496
497	# # TODO: Plane Elasticity
498	# if Set.InPlaneElasticity:
499	# gt, Kt, EBulk = KgBulk(geo_0, Geo, Set)
500	# K += Kt
501	# g += gt
502	# E += EBulk
503	# Energies["Bulk"] = EBulk
504
505	# Bending Energy
506	# TODO
507
508	# Triangle Energy Barrier
509	if Set.EnergyBarrierA:	×
510	kg_Tri = KgTriEnergyBarrier(Geo)	×
511	kg_Tri.compute_work(Geo, Set, None, False)	×
512	g += kg_Tri.g	×
513	energies["TriEnergyBarrier"] = kg_Tri.energy	×
514
515	# Triangle Energy Barrier Aspect Ratio
516	if Set.EnergyBarrierAR:	×
517	kg_TriAR = KgTriAREnergyBarrier(Geo)	×
518	kg_TriAR.compute_work(Geo, Set, None, False)	×
519	g += kg_TriAR.g	×
520	energies["TriEnergyBarrierAR"] = kg_TriAR.energy	×
521
522	# Propulsion Forces
523	# TODO
524
525	# Contractility
526	if Set.Contractility:	×
527	kg_lt = KgContractility(Geo)	×
528	kg_lt.compute_work(Geo, Set, None, False)	×
529	g += kg_lt.g	×
530	energies["Contractility"] = kg_lt.energy	×
531
532	# Contractility as external force
533	if Set.Contractility_external:	×
534	kg_ext = KgContractilityExternal(Geo)	×
535	kg_ext.compute_work(Geo, Set, None, False)	×
536	g += kg_ext.g	×
537	energies["Contractility_external"] = kg_ext.energy	×
538
539	# Substrate
540	if Set.Substrate == 2:	×
541	kg_subs = KgSubstrate(Geo)	×
542	kg_subs.compute_work(Geo, Set, None, False)	×
543	g += kg_subs.g	×
544	energies["Substrate"] = kg_subs.energy	×
545
546	return g, energies	×
547
548
UNCOV 549	def remeshing_cells(Geo_0, Geo_n, Geo, Dofs, Set, cells_to_change, ghost_node):	×
550	"""
551	Newton-Raphson method
552	:param Geo_0:
553	:param Geo_n:
554	:param Geo:
555	:param Dofs:
556	:param Set:s
557	:return:
558	"""
559
560	def objective_function(dy, dof):	×
561	# Reshape dy to match the geometry of the system
562	dy_reshaped = np.reshape(dy, (Geo.numF + Geo.numY + Geo.nCells, 3))	×
563
564	# Create a copy of Geo to not change the original one
565	Geo_copy = Geo.copy()	×
566
567	# Update the vertices in the copy
568	Geo_copy.update_vertices(dy_reshaped)	×
569
570	# Update the measures in the copy
571	Geo_copy.update_measures()	×
572
573	# Compute the energies
574	g, _, energy_total, _ = KgGlobal(Geo_0, Geo_n, Geo_copy, Set)	×
575
576	logger.info('Energy: ' + str(energy_total))	×
577	gr = np.linalg.norm(g[dof])	×
578	logger.info('\|\|gr\|\|= ' + str(gr))	×
579
580	# Return the total energy
581	return gr	×
582
583	def gradient_function(dy, dof):	×
584	# Reshape dy to match the geometry of the system
585	dy_reshaped = np.reshape(dy, (Geo.numF + Geo.numY + Geo.nCells, 3))	×
586
587	# Create a copy of Geo to not change the original one
588	Geo_copy = Geo.copy()	×
589
590	# Update the vertices in the copy
591	Geo_copy.update_vertices(dy_reshaped)	×
592
593	# Update the measures in the copy
594	Geo_copy.update_measures()	×
595
596	# Compute the energies
597	g, _, _, _ = KgGlobal(Geo_0, Geo_n, Geo_copy, Set)	×
598
599	g[dof == False] = 0	×
600
601	# Return the gradient
602	return g.flatten()	×
603
604	Geo.remeshing = True	×
605
606	if np.isin(ghost_node, Geo.XgBottom).any():	×
607	interface_type = 2	×
608	elif np.isin(ghost_node, Geo.XgTop).any():	×
609	interface_type = 0	×
610
611	Dofs.get_remeshing_dofs(Geo, cells_to_change, interface_type)	×
612
613	# Set the degrees of freedom
614	dof = Dofs.remeshing	×
615
616	dy_initial = np.zeros((Geo.numY + Geo.numF + Geo.nCells) * 3, dtype=np.float64)	×
617	g, K, _, _ = KgGlobal(Geo_0, Geo_n, Geo, Set)	×
618	dy_initial[dof] = ml_divide(K, dof, g)	×
619
620	lower_bounds = np.zeros_like(dy_initial)	×
621	upper_bounds = 10e-1 * np.ones_like(dy_initial)	×
622	bounds = Bounds(lower_bounds, upper_bounds)	×
623
624	result = minimize(	×
625	fun=objective_function, # Defined as before
626	x0=dy_initial,
627	args=(dof,),
628	method='L-BFGS-B', # Newton-CG
629	jac=gradient_function,
630	bounds=bounds,
631	options={'disp': True}
632	)
633
634	dy_reshaped = np.reshape(result.x, (Geo.numF + Geo.numY + Geo.nCells, 3))	×
635	return dy_reshaped	×

Pablo1990 / pyVertexModel / 11363842044

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous