11234788583

Committed 08 Oct 2024 11:36AM UTC coverage: 8.778% (-0.004%) from 8.782%

Build # 11234788583

Build Type

push

github

Committed by

Pablo1990

Commit Message

it is not perfect, but it is working

Run Details

0 of 1007 branches covered (0.0%)

Branch coverage included in aggregate %.

0 of 3 new or added lines in 1 file covered. (0.0%)

1 existing line in 1 file now uncovered.

563 of 5407 relevant lines covered (10.41%)

0.1 hits per line

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

6.07

/src/pyVertexModel/algorithm/newtonRaphson.py

import logging

import numpy as np
from scipy.optimize import minimize, Bounds

from src.pyVertexModel.Kg.kgContractility import KgContractility
from src.pyVertexModel.Kg.kgContractility_external import KgContractilityExternal
from src.pyVertexModel.Kg.kgSubstrate import KgSubstrate
from src.pyVertexModel.Kg.kgSurfaceCellBasedAdhesion import KgSurfaceCellBasedAdhesion
from src.pyVertexModel.Kg.kgTriAREnergyBarrier import KgTriAREnergyBarrier
from src.pyVertexModel.Kg.kgTriEnergyBarrier import KgTriEnergyBarrier
from src.pyVertexModel.Kg.kgViscosity import KgViscosity
from src.pyVertexModel.Kg.kgVolume import KgVolume

logger = logging.getLogger("pyVertexModel")


def solve_remodeling_step(geo_0, geo_n, geo, dofs, c_set):
    """
    This function solves local problem to obtain the position of the newly remodeled vertices with prescribed settings
    (Set.***_LP), e.g. Set.lambda_LP.
    :param geo_0:
    :param geo_n:
    :param geo:
    :param dofs:
    :param c_set:
    :return:
    """

    logger.info('=====>> Solving Local Problem....')
    geo.remodelling = True
    original_nu = c_set.nu
    original_nu0 = c_set.nu0
    original_lambdaB = c_set.lambdaB
    original_lambdaR = c_set.lambdaR

    c_set.MaxIter = c_set.MaxIter0 * 3
    c_set.lambdaB = original_lambdaB * 2
    c_set.lambdaR = original_lambdaR * 0.0001

    for n_id, nu_factor in enumerate(np.linspace(10, 1, 3)):
        c_set.nu0 = original_nu0 * nu_factor
        c_set.nu = original_nu * nu_factor
        g, k, _, _ = KgGlobal(geo_0, geo_n, geo, c_set)

        dy = np.zeros(((geo.numY + geo.numF + geo.nCells) * 3, 1), dtype=np.float64)
        dyr = np.linalg.norm(dy[dofs.remodel, 0])
        gr = np.linalg.norm(g[dofs.remodel])
        logger.info(
            f'Local Problem ->Iter: 0, ||gr||= {gr:.3e} ||dyr||= {dyr:.3e}  nu/nu0={c_set.nu / c_set.nu0:.3e}  '
            f'dt/dt0={c_set.dt / c_set.dt0:.3g}')

        geo, g, k, energy, c_set, gr, dyr, dy = newton_raphson(geo_0, geo_n, geo, dofs, c_set, k, g, -1, -1)

        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])):
            logger.info(f'Local Problem did not converge after {c_set.iter} iterations.')
            has_converged = False
            break
        else:
            logger.info(f'=====>> Local Problem converged in {c_set.iter} iterations.')
            has_converged = True

    geo.remodelling = False

    c_set.MaxIter = c_set.MaxIter0
    c_set.nu = original_nu
    c_set.nu0 = original_nu0
    c_set.lambdaB = original_lambdaB
    c_set.lambdaR = original_lambdaR

    return geo, c_set, has_converged


def newton_raphson(Geo_0, Geo_n, Geo, Dofs, Set, K, g, numStep, t, implicit_method=True):
    """
    Newton-Raphson method
    :param Geo_0:
    :param Geo_n:
    :param Geo:
    :param Dofs:
    :param Set:
    :param K:
    :param g:
    :param numStep:
    :param t:
    :return:
    """
    if Geo.remodelling:
        dof = Dofs.remodel
    else:
        dof = Dofs.Free

    dy = np.zeros(((Geo.numY + Geo.numF + Geo.nCells) * 3, 1), dtype=np.float64)
    dyr = np.linalg.norm(dy[dof, 0])
    gr = np.linalg.norm(g[dof])
    gr0 = gr

    logger.info(f"Step: {numStep}, Iter: 0 ||gr||= {gr} ||dyr||= {dyr} dt/dt0={Set.dt / Set.dt0:.3g}")

    Set.iter = 1
    auxgr = np.zeros(3, dtype=np.float64)
    auxgr[0] = gr
    ig = 0
    Energy = 0

    if implicit_method is True:
        while (gr > Set.tol or dyr > Set.tol) and Set.iter < Set.MaxIter:
            Energy, K, dyr, g, gr, ig, auxgr, dy = newton_raphson_iteration(Dofs, Geo, Geo_0, Geo_n, K, Set, auxgr, dof,
                                                                            dy, g, gr0, ig, numStep, t)
    else:
        Geo, dy, gr = newton_raphson_iteration_explicit(Geo, Set, dof, dy, g)
        dyr = 0

    return Geo, g, K, Energy, Set, gr, dyr, dy


def newton_raphson_iteration(Dofs, Geo, Geo_0, Geo_n, K, Set, aux_gr, dof, dy, g, gr0, ig, numStep, t):
    """
    Perform a single iteration of the Newton-Raphson method for solving a system of nonlinear equations.

    Parameters:
    Dofs (object): Object containing information about the degrees of freedom in the system.
    Geo (object): Object containing the current geometry of the system.
    Geo_0 (object): Object containing the initial geometry of the system.
    Geo_n (object): Object containing the geometry of the system at the previous time step.
    K (ndarray): The Jacobian matrix of the system.
    Set (object): Object containing various settings for the simulation.
    aux_gr (ndarray): Array containing the norms of the gradient at the previous three steps.
    dof (ndarray): Array containing the indices of the free degrees of freedom in the system.
    dy (ndarray): Array containing the current guess for the solution of the system.
    g (ndarray): The gradient of the system.
    gr0 (float): The norm of the gradient at the start of the current time step.
    ig (int): Index used to cycle through the elements of aux_gr.
    numStep (int): The current time step number.
    t (float): The current time.

    Returns:
    energy_total (float): The total energy of the system after the current iteration.
    K (ndarray): The updated Jacobian matrix of the system.
    dyr (float): The norm of the change in the solution guess during the current iteration.
    g (ndarray): The updated gradient of the system.
    gr (float): The norm of the updated gradient of the system.
    ig (int): The updated index used to cycle through the elements of aux_gr.
    aux_gr (ndarray): The updated array containing the norms of the gradient at the previous three steps.
    dy (ndarray): The updated guess for the solution of the system.
    """
    dy[dof, 0] = ml_divide(K, dof, g)

    alpha = line_search(Geo_0, Geo_n, Geo, dof, Set, g, dy)
    dy_reshaped = np.reshape(dy * alpha, (Geo.numF + Geo.numY + Geo.nCells, 3))
    Geo.update_vertices(dy_reshaped)
    Geo.update_measures()
    g, K, energy_total, _ = KgGlobal(Geo_0, Geo_n, Geo, Set)

    dyr = np.linalg.norm(dy[dof, 0])
    gr = np.linalg.norm(g[dof])
    logger.info(f"Step: {numStep}, Iter: {Set.iter}, Time: {t} ||gr||= {gr:.3e} ||dyr||= {dyr:.3e} alpha= {alpha:.3e}"
                f" nu/nu0={Set.nu / Set.nu0:.3g}")
    Set.iter += 1

    # Checking if the three previous steps are very similar. Thus, the solution is not converging
    aux_gr[ig] = gr
    ig = 0 if ig == 2 else ig + 1
    if (
            all(aux_gr[i] != 0 for i in range(3)) and
            abs(aux_gr[0] - aux_gr[1]) / aux_gr[0] < 1e-3 and
            abs(aux_gr[0] - aux_gr[2]) / aux_gr[0] < 1e-3 and
            abs(aux_gr[2] - aux_gr[1]) / aux_gr[2] < 1e-3
    ) or (
            gr0 != 0 and abs((gr0 - gr) / gr0) > 1e3
    ):
        Set.iter = Set.MaxIter

    return energy_total, K, dyr, g, gr, ig, aux_gr, dy


def newton_raphson_iteration_explicit(Geo, Set, dof, dy, g):
    """
    Explicit update method
    :param Geo_0:
    :param Geo_n:
    :param Geo:
    :param Dofs:
    :param Set:
    :return:
    """
    dy[dof, 0] = -Set.dt / Set.nu * g[dof]

    # Update border ghost nodes with the same displacement as the corresponding real nodes
    for numCell in Geo.BorderCells:
        cell = Geo.Cells[numCell]

        tets_to_check = cell.T[np.any(np.isin(cell.T, Geo.BorderGhostNodes), axis=1)]
        global_ids_to_update = cell.globalIds[np.any(np.isin(cell.T, Geo.BorderGhostNodes), axis=1)]

        for i, tet in enumerate(tets_to_check):
            for j, node in enumerate(tet):
                if node in Geo.BorderGhostNodes:
                    global_id = global_ids_to_update[i]
                    # Iterate over each element in tet and access the opposite_cell attribute
                    opposite_cell = Geo.Cells[Geo.Cells[int(node)].opposite_cell]
                    opposite_cells = np.unique([Geo.Cells[int(t)].opposite_cell for t in tet if Geo.Cells[int(t)].opposite_cell is not None])

                    # Get the most similar tetrahedron in the opposite cell
                    if len(opposite_cells) < 3:
                        possible_tets = opposite_cell.T[[np.sum(np.isin(opposite_cells, tet)) == len(opposite_cells) for tet in opposite_cell.T]]
                    else:
                        possible_tets = opposite_cell.T[[np.sum(np.isin(opposite_cells, tet)) >= 2 for tet in opposite_cell.T]]

                    if possible_tets.shape[0] == 0:
                        possible_tets = opposite_cell.T[[np.sum(np.isin(opposite_cells, tet)) == 1 for tet in opposite_cell.T]]

                    if possible_tets.shape[0] > 0:
                        # Filter the possible tets to get its correct location (XgTop, XgBottom or XgLateral)
                        if possible_tets.shape[0] > 1:
                            old_possible_tets = possible_tets
                            scutoid = False
                            if np.isin(tet, Geo.XgTop).any():
                                possible_tets = possible_tets[np.isin(possible_tets, Geo.XgTop).any(axis=1)]
                            elif np.isin(tet, Geo.XgBottom).any():
                                possible_tets = possible_tets[np.isin(possible_tets, Geo.XgBottom).any(axis=1)]
                            else:
                                # Get the tets that are not in the top or bottom
                                possible_tets = possible_tets[np.logical_not(np.isin(possible_tets, Geo.XgTop).any(axis=1))]
                                scutoid = True

                            if possible_tets.shape[0] == 0:
                                possible_tets = old_possible_tets

                        # Compare tets with the number of ghost nodes
                        if possible_tets.shape[0] > 1 and not scutoid:
                            old_possible_tets = possible_tets
                            # Get the tets that have the same number of ghost nodes as the original tet
                            possible_tets = possible_tets[np.sum(np.isin(tet[tet!=node], Geo.XgID)) == np.sum(np.isin(possible_tets, Geo.XgID), axis=1)]

                            if possible_tets.shape[0] == 0:
                                possible_tets = old_possible_tets

                        if possible_tets.shape[0] > 1 and not scutoid:
                            old_possible_tets = possible_tets
                            # Get the tet that has the same number of ghost nodes as the original tet
                            possible_tets = possible_tets[np.sum(np.isin(tet[tet!=node], Geo.XgID)) == np.sum(np.isin(possible_tets, np.concatenate([Geo.XgTop, Geo.XgBottom])), axis=1)]

                            if possible_tets.shape[0] == 0:
                                possible_tets = old_possible_tets
                    else:
                        print('Error in Tet ID: ', tet)

                    if possible_tets.shape[0] == 1:
                        opposite_global_id = opposite_cell.globalIds[np.where(np.all(opposite_cell.T == possible_tets[0], axis=1))[0][0]]
                        dy[global_id * 3:global_id * 3 + 3, 0] = dy[opposite_global_id * 3:opposite_global_id * 3 + 3, 0]
                    elif possible_tets.shape[0] > 1:
                        # TODO: It is either the first one or an average of the two
                        if scutoid:
                            avg_dy = np.zeros(3)
                            # Average of the dys
                            for c_tet in possible_tets:
                                opposite_global_id = opposite_cell.globalIds[
                                    np.where(np.all(opposite_cell.T == c_tet, axis=1))[0][0]]
                                avg_dy += dy[opposite_global_id * 3:opposite_global_id * 3 + 3, 0]
                            avg_dy /= possible_tets.shape[0]
                            dy[global_id * 3:global_id * 3 + 3, 0] = avg_dy
                        else:
                            opposite_global_id = opposite_cell.globalIds[
                                np.where(np.all(opposite_cell.T == possible_tets[0], axis=1))[0][0]]
                            dy[global_id * 3:global_id * 3 + 3, 0] = dy[opposite_global_id * 3:opposite_global_id * 3 + 3, 0]

    dy_reshaped = np.reshape(dy, (Geo.numF + Geo.numY + Geo.nCells, 3))
    Geo.update_vertices(dy_reshaped)
    Geo.update_measures()

    g, energies = gGlobal(Geo, Geo, Geo, Set, Set.implicit_method)
    gr = np.linalg.norm(g[dof])
    Set.iter = Set.MaxIter

    return Geo, dy, gr


def ml_divide(K, dof, g):
    """
    Solve the linear system K * dy = g
    :param K:
    :param dof:
    :param g:
    :return:
    """
    # dy[dof] = kg_functions.mldivide_np(K[np.ix_(dof, dof)], g[dof])
    return -np.linalg.solve(K[np.ix_(dof, dof)], g[dof])


def line_search(Geo_0, Geo_n, geo, dof, Set, gc, dy):
    """
    Line search to find the best alpha to minimize the energy
    :param Geo_0:   Initial geometry.
    :param Geo_n:   Geometry at the previous time step
    :param geo:     Current geometry
    :param Dofs:    Dofs object
    :param Set:     Set object
    :param gc:      Gradient at the current step
    :param dy:      Displacement at the current step
    :return:        alpha
    """
    dy_reshaped = np.reshape(dy, (geo.numF + geo.numY + geo.nCells, 3))

    # Create a copy of geo to not change the original one
    Geo_copy = geo.copy()

    Geo_copy.update_vertices(dy_reshaped)
    Geo_copy.update_measures()

    g, _ = gGlobal(Geo_0, Geo_n, Geo_copy, Set)

    gr0 = np.linalg.norm(gc[dof])
    gr = np.linalg.norm(g[dof])

    if gr0 < gr:
        R0 = np.dot(dy[dof, 0], gc[dof])
        R1 = np.dot(dy[dof, 0], g[dof])

        R = R0 / R1
        alpha1 = (R / 2) + np.sqrt((R / 2) ** 2 - R)
        alpha2 = (R / 2) - np.sqrt((R / 2) ** 2 - R)

        if np.isreal(alpha1) and 2 > alpha1 > 1e-3:
            alpha = alpha1
        elif np.isreal(alpha2) and 2 > alpha2 > 1e-3:
            alpha = alpha2
        else:
            alpha = 0.1
    else:
        alpha = 1

    return alpha


def KgGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True):
    """
    Compute the global Jacobian matrix and the global gradient
    :param Geo_0:
    :param Geo_n:
    :param Geo:
    :param Set:
    :param implicit_method:
    :return:
    """

    # Surface Energy
    kg_SA = KgSurfaceCellBasedAdhesion(Geo)
    kg_SA.compute_work(Geo, Set)
    g = kg_SA.g
    K = kg_SA.K
    energy_total = kg_SA.energy
    energies = {"Surface": kg_SA.energy}

    # Volume Energy
    kg_Vol = KgVolume(Geo)
    kg_Vol.compute_work(Geo, Set)
    K += kg_Vol.K
    g += kg_Vol.g
    energy_total += kg_Vol.energy
    energies["Volume"] = kg_Vol.energy

    # # TODO: Plane Elasticity
    # if Set.InPlaneElasticity:
    #     gt, Kt, EBulk = KgBulk(geo_0, Geo, Set)
    #     K += Kt
    #     g += gt
    #     E += EBulk
    #     Energies["Bulk"] = EBulk

    # Bending Energy
    # TODO

    # Propulsion Forces
    # TODO

    # Contractility
    if Set.Contractility:
        kg_lt = KgContractility(Geo)
        kg_lt.compute_work(Geo, Set)
        g += kg_lt.g
        K += kg_lt.K
        energy_total += kg_lt.energy
        energies["Contractility"] = kg_lt.energy

    if Set.Contractility_external:
        kg_ext = KgContractilityExternal(Geo)
        kg_ext.compute_work(Geo, Set)
        g += kg_ext.g
        K += kg_ext.K
        energy_total += kg_ext.energy
        energies["Contractility_external"] = kg_ext.energy

    # Substrate
    if Set.Substrate == 2:
        kg_subs = KgSubstrate(Geo)
        kg_subs.compute_work(Geo, Set)
        g += kg_subs.g
        K += kg_subs.K
        energy_total += kg_subs.energy
        energies["Substrate"] = kg_subs.energy

    # Triangle Energy Barrier
    if Set.EnergyBarrierA:
        kg_Tri = KgTriEnergyBarrier(Geo)
        kg_Tri.compute_work(Geo, Set)
        g += kg_Tri.g
        K += kg_Tri.K
        energy_total += kg_Tri.energy
        energies["TriEnergyBarrier"] = kg_Tri.energy

    # Triangle Energy Barrier Aspect Ratio
    if Set.EnergyBarrierAR:
        kg_TriAR = KgTriAREnergyBarrier(Geo)
        kg_TriAR.compute_work(Geo, Set)
        g += kg_TriAR.g
        K += kg_TriAR.K
        energy_total += kg_TriAR.energy
        energies["TriEnergyBarrierAR"] = kg_TriAR.energy

    if implicit_method is True:
        # Viscous Energy
        kg_Viscosity = KgViscosity(Geo)
        kg_Viscosity.compute_work(Geo, Set, Geo_n)
        g += kg_Viscosity.g
        K += kg_Viscosity.K
        energy_total = kg_Viscosity.energy
        energies["Viscosity"] = kg_Viscosity.energy

    return g, K, energy_total, energies


def gGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True):
    """
    Compute the global gradient
    :param Geo_0:
    :param Geo_n:
    :param Geo:
    :param Set:
    :param implicit_method:
    :return:
    """
    g = np.zeros((Geo.numY + Geo.numF + Geo.nCells) * 3, dtype=np.float64)
    energies = {}

    # Surface Energy
    if Set.lambdaS1 > 0:
        kg_SA = KgSurfaceCellBasedAdhesion(Geo)
        kg_SA.compute_work(Geo, Set, None, False)
        g += kg_SA.g
        energies["Surface"] = kg_SA.energy

    # Volume Energy
    if Set.lambdaV > 0:
        kg_Vol = KgVolume(Geo)
        kg_Vol.compute_work(Geo, Set, None, False)
        g += kg_Vol.g[:]
        energies["Volume"] = kg_Vol.energy

    if implicit_method is True:
        # Viscous Energy
        kg_Viscosity = KgViscosity(Geo)
        kg_Viscosity.compute_work(Geo, Set, Geo_n, False)
        g += kg_Viscosity.g
        energies["Viscosity"] = kg_Viscosity.energy

    # # TODO: Plane Elasticity
    # if Set.InPlaneElasticity:
    #     gt, Kt, EBulk = KgBulk(geo_0, Geo, Set)
    #     K += Kt
    #     g += gt
    #     E += EBulk
    #     Energies["Bulk"] = EBulk

    # Bending Energy
    # TODO

    # Triangle Energy Barrier
    if Set.EnergyBarrierA:
        kg_Tri = KgTriEnergyBarrier(Geo)
        kg_Tri.compute_work(Geo, Set, None, False)
        g += kg_Tri.g
        energies["TriEnergyBarrier"] = kg_Tri.energy

    # Triangle Energy Barrier Aspect Ratio
    if Set.EnergyBarrierAR:
        kg_TriAR = KgTriAREnergyBarrier(Geo)
        kg_TriAR.compute_work(Geo, Set, None, False)
        g += kg_TriAR.g
        energies["TriEnergyBarrierAR"] = kg_TriAR.energy

    # Propulsion Forces
    # TODO

    # Contractility
    if Set.Contractility:
        kg_lt = KgContractility(Geo)
        kg_lt.compute_work(Geo, Set, None, False)
        g += kg_lt.g
        energies["Contractility"] = kg_lt.energy

    # Contractility as external force
    if Set.Contractility_external:
        kg_ext = KgContractilityExternal(Geo)
        kg_ext.compute_work(Geo, Set, None, False)
        g += kg_ext.g
        energies["Contractility_external"] = kg_ext.energy

    # Substrate
    if Set.Substrate == 2:
        kg_subs = KgSubstrate(Geo)
        kg_subs.compute_work(Geo, Set, None, False)
        g += kg_subs.g
        energies["Substrate"] = kg_subs.energy

    return g, energies


def remeshing_cells(Geo_0, Geo_n, Geo, Dofs, Set, cells_to_change, ghost_node):
    """
    Newton-Raphson method
    :param Geo_0:
    :param Geo_n:
    :param Geo:
    :param Dofs:
    :param Set:s
    :return:
    """

    def objective_function(dy, dof):
        # Reshape dy to match the geometry of the system
        dy_reshaped = np.reshape(dy, (Geo.numF + Geo.numY + Geo.nCells, 3))

        # Create a copy of Geo to not change the original one
        Geo_copy = Geo.copy()

        # Update the vertices in the copy
        Geo_copy.update_vertices(dy_reshaped)

        # Update the measures in the copy
        Geo_copy.update_measures()

        # Compute the energies
        g, _, energy_total, _ = KgGlobal(Geo_0, Geo_n, Geo_copy, Set)

        logger.info('Energy: ' + str(energy_total))
        gr = np.linalg.norm(g[dof])
        logger.info('||gr||= ' + str(gr))

        # Return the total energy
        return gr

    def gradient_function(dy, dof):
        # Reshape dy to match the geometry of the system
        dy_reshaped = np.reshape(dy, (Geo.numF + Geo.numY + Geo.nCells, 3))

        # Create a copy of Geo to not change the original one
        Geo_copy = Geo.copy()

        # Update the vertices in the copy
        Geo_copy.update_vertices(dy_reshaped)

        # Update the measures in the copy
        Geo_copy.update_measures()

        # Compute the energies
        g, _, _, _ = KgGlobal(Geo_0, Geo_n, Geo_copy, Set)

        g[dof == False] = 0

        # Return the gradient
        return g.flatten()

    Geo.remeshing = True

    if np.isin(ghost_node, Geo.XgBottom).any():
        interface_type = 2
    elif np.isin(ghost_node, Geo.XgTop).any():
        interface_type = 0

    Dofs.get_remeshing_dofs(Geo, cells_to_change, interface_type)

    # Set the degrees of freedom
    dof = Dofs.remeshing

    dy_initial = np.zeros((Geo.numY + Geo.numF + Geo.nCells) * 3, dtype=np.float64)
    g, K, _, _ = KgGlobal(Geo_0, Geo_n, Geo, Set)
    dy_initial[dof] = ml_divide(K, dof, g)

    lower_bounds = np.zeros_like(dy_initial)
    upper_bounds = 10e-1 * np.ones_like(dy_initial)
    bounds = Bounds(lower_bounds, upper_bounds)

    result = minimize(
        fun=objective_function,  # Defined as before
        x0=dy_initial,
        args=(dof,),
        method='L-BFGS-B',  # Newton-CG
        jac=gradient_function,
        bounds=bounds,
        options={'disp': True}
    )

    dy_reshaped = np.reshape(result.x, (Geo.numF + Geo.numY + Geo.nCells, 3))
    return dy_reshaped

1	import logging	1✔
2
3	import numpy as np	1✔
4	from scipy.optimize import minimize, Bounds	1✔
5
6	from src.pyVertexModel.Kg.kgContractility import KgContractility	1✔
7	from src.pyVertexModel.Kg.kgContractility_external import KgContractilityExternal	1✔
8	from src.pyVertexModel.Kg.kgSubstrate import KgSubstrate	1✔
9	from src.pyVertexModel.Kg.kgSurfaceCellBasedAdhesion import KgSurfaceCellBasedAdhesion	1✔
10	from src.pyVertexModel.Kg.kgTriAREnergyBarrier import KgTriAREnergyBarrier	1✔
11	from src.pyVertexModel.Kg.kgTriEnergyBarrier import KgTriEnergyBarrier	1✔
12	from src.pyVertexModel.Kg.kgViscosity import KgViscosity	1✔
13	from src.pyVertexModel.Kg.kgVolume import KgVolume	1✔
14
15	logger = logging.getLogger("pyVertexModel")	1✔
16
17
18	def solve_remodeling_step(geo_0, geo_n, geo, dofs, c_set):	1✔
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
74	def newton_raphson(Geo_0, Geo_n, Geo, Dofs, Set, K, g, numStep, t, implicit_method=True):	1✔
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
117	def newton_raphson_iteration(Dofs, Geo, Geo_0, Geo_n, K, Set, aux_gr, dof, dy, g, gr0, ig, numStep, t):	1✔
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
177	def newton_raphson_iteration_explicit(Geo, Set, dof, dy, g):	1✔
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	for numCell in Geo.BorderCells:	×
191	cell = Geo.Cells[numCell]	×
192
193	tets_to_check = cell.T[np.any(np.isin(cell.T, Geo.BorderGhostNodes), axis=1)]	×
194	global_ids_to_update = cell.globalIds[np.any(np.isin(cell.T, Geo.BorderGhostNodes), axis=1)]	×
195
196	for i, tet in enumerate(tets_to_check):	×
197	for j, node in enumerate(tet):	×
198	if node in Geo.BorderGhostNodes:	×
199	global_id = global_ids_to_update[i]	×
200	# Iterate over each element in tet and access the opposite_cell attribute
201	opposite_cell = Geo.Cells[Geo.Cells[int(node)].opposite_cell]	×
202	opposite_cells = np.unique([Geo.Cells[int(t)].opposite_cell for t in tet if Geo.Cells[int(t)].opposite_cell is not None])	×
203
204	# Get the most similar tetrahedron in the opposite cell
205	if len(opposite_cells) < 3:	×
206	possible_tets = opposite_cell.T[[np.sum(np.isin(opposite_cells, tet)) == len(opposite_cells) for tet in opposite_cell.T]]	×
207	else:
208	possible_tets = opposite_cell.T[[np.sum(np.isin(opposite_cells, tet)) >= 2 for tet in opposite_cell.T]]	×
209
NEW 210	if possible_tets.shape[0] == 0:	×
NEW 211	possible_tets = opposite_cell.T[[np.sum(np.isin(opposite_cells, tet)) == 1 for tet in opposite_cell.T]]	×
212
UNCOV 213	if possible_tets.shape[0] > 0:	×
214	# Filter the possible tets to get its correct location (XgTop, XgBottom or XgLateral)
215	if possible_tets.shape[0] > 1:	×
216	old_possible_tets = possible_tets	×
217	scutoid = False	×
218	if np.isin(tet, Geo.XgTop).any():	×
219	possible_tets = possible_tets[np.isin(possible_tets, Geo.XgTop).any(axis=1)]	×
220	elif np.isin(tet, Geo.XgBottom).any():	×
221	possible_tets = possible_tets[np.isin(possible_tets, Geo.XgBottom).any(axis=1)]	×
222	else:
223	# Get the tets that are not in the top or bottom
224	possible_tets = possible_tets[np.logical_not(np.isin(possible_tets, Geo.XgTop).any(axis=1))]	×
225	scutoid = True	×
226
227	if possible_tets.shape[0] == 0:	×
228	possible_tets = old_possible_tets	×
229
230	# Compare tets with the number of ghost nodes
231	if possible_tets.shape[0] > 1 and not scutoid:	×
232	old_possible_tets = possible_tets	×
233	# Get the tets that have the same number of ghost nodes as the original tet
234	possible_tets = possible_tets[np.sum(np.isin(tet[tet!=node], Geo.XgID)) == np.sum(np.isin(possible_tets, Geo.XgID), axis=1)]	×
235
236	if possible_tets.shape[0] == 0:	×
237	possible_tets = old_possible_tets	×
238
239	if possible_tets.shape[0] > 1 and not scutoid:	×
240	old_possible_tets = possible_tets	×
241	# Get the tet that has the same number of ghost nodes as the original tet
242	possible_tets = possible_tets[np.sum(np.isin(tet[tet!=node], Geo.XgID)) == np.sum(np.isin(possible_tets, np.concatenate([Geo.XgTop, Geo.XgBottom])), axis=1)]	×
243
244	if possible_tets.shape[0] == 0:	×
245	possible_tets = old_possible_tets	×
246	else:
NEW 247	print('Error in Tet ID: ', tet)	×
248
249	if possible_tets.shape[0] == 1:	×
250	opposite_global_id = opposite_cell.globalIds[np.where(np.all(opposite_cell.T == possible_tets[0], axis=1))[0][0]]	×
251	dy[global_id * 3:global_id * 3 + 3, 0] = dy[opposite_global_id * 3:opposite_global_id * 3 + 3, 0]	×
252	elif possible_tets.shape[0] > 1:	×
253	# TODO: It is either the first one or an average of the two
254	if scutoid:	×
255	avg_dy = np.zeros(3)	×
256	# Average of the dys
257	for c_tet in possible_tets:	×
258	opposite_global_id = opposite_cell.globalIds[	×
259	np.where(np.all(opposite_cell.T == c_tet, axis=1))[0][0]]
260	avg_dy += dy[opposite_global_id * 3:opposite_global_id * 3 + 3, 0]	×
261	avg_dy /= possible_tets.shape[0]	×
262	dy[global_id * 3:global_id * 3 + 3, 0] = avg_dy	×
263	else:
264	opposite_global_id = opposite_cell.globalIds[	×
265	np.where(np.all(opposite_cell.T == possible_tets[0], axis=1))[0][0]]
266	dy[global_id * 3:global_id * 3 + 3, 0] = dy[opposite_global_id * 3:opposite_global_id * 3 + 3, 0]	×
267
268	dy_reshaped = np.reshape(dy, (Geo.numF + Geo.numY + Geo.nCells, 3))	×
269	Geo.update_vertices(dy_reshaped)	×
270	Geo.update_measures()	×
271
272	g, energies = gGlobal(Geo, Geo, Geo, Set, Set.implicit_method)	×
273	gr = np.linalg.norm(g[dof])	×
274	Set.iter = Set.MaxIter	×
275
276	return Geo, dy, gr	×
277
278
279	def ml_divide(K, dof, g):	1✔
280	"""
281	Solve the linear system K * dy = g
282	:param K:
283	:param dof:
284	:param g:
285	:return:
286	"""
287	# dy[dof] = kg_functions.mldivide_np(K[np.ix_(dof, dof)], g[dof])
288	return -np.linalg.solve(K[np.ix_(dof, dof)], g[dof])	×
289
290
291	def line_search(Geo_0, Geo_n, geo, dof, Set, gc, dy):	1✔
292	"""
293	Line search to find the best alpha to minimize the energy
294	:param Geo_0: Initial geometry.
295	:param Geo_n: Geometry at the previous time step
296	:param geo: Current geometry
297	:param Dofs: Dofs object
298	:param Set: Set object
299	:param gc: Gradient at the current step
300	:param dy: Displacement at the current step
301	:return: alpha
302	"""
303	dy_reshaped = np.reshape(dy, (geo.numF + geo.numY + geo.nCells, 3))	×
304
305	# Create a copy of geo to not change the original one
306	Geo_copy = geo.copy()	×
307
308	Geo_copy.update_vertices(dy_reshaped)	×
309	Geo_copy.update_measures()	×
310
311	g, _ = gGlobal(Geo_0, Geo_n, Geo_copy, Set)	×
312
313	gr0 = np.linalg.norm(gc[dof])	×
314	gr = np.linalg.norm(g[dof])	×
315
316	if gr0 < gr:	×
317	R0 = np.dot(dy[dof, 0], gc[dof])	×
318	R1 = np.dot(dy[dof, 0], g[dof])	×
319
320	R = R0 / R1	×
321	alpha1 = (R / 2) + np.sqrt((R / 2) ** 2 - R)	×
322	alpha2 = (R / 2) - np.sqrt((R / 2) ** 2 - R)	×
323
324	if np.isreal(alpha1) and 2 > alpha1 > 1e-3:	×
325	alpha = alpha1	×
326	elif np.isreal(alpha2) and 2 > alpha2 > 1e-3:	×
327	alpha = alpha2	×
328	else:
329	alpha = 0.1	×
330	else:
331	alpha = 1	×
332
333	return alpha	×
334
335
336	def KgGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True):	1✔
337	"""
338	Compute the global Jacobian matrix and the global gradient
339	:param Geo_0:
340	:param Geo_n:
341	:param Geo:
342	:param Set:
343	:param implicit_method:
344	:return:
345	"""
346
347	# Surface Energy
348	kg_SA = KgSurfaceCellBasedAdhesion(Geo)	×
349	kg_SA.compute_work(Geo, Set)	×
350	g = kg_SA.g	×
351	K = kg_SA.K	×
352	energy_total = kg_SA.energy	×
353	energies = {"Surface": kg_SA.energy}	×
354
355	# Volume Energy
356	kg_Vol = KgVolume(Geo)	×
357	kg_Vol.compute_work(Geo, Set)	×
358	K += kg_Vol.K	×
359	g += kg_Vol.g	×
360	energy_total += kg_Vol.energy	×
361	energies["Volume"] = kg_Vol.energy	×
362
363	# # TODO: Plane Elasticity
364	# if Set.InPlaneElasticity:
365	# gt, Kt, EBulk = KgBulk(geo_0, Geo, Set)
366	# K += Kt
367	# g += gt
368	# E += EBulk
369	# Energies["Bulk"] = EBulk
370
371	# Bending Energy
372	# TODO
373
374	# Propulsion Forces
375	# TODO
376
377	# Contractility
378	if Set.Contractility:	×
379	kg_lt = KgContractility(Geo)	×
380	kg_lt.compute_work(Geo, Set)	×
381	g += kg_lt.g	×
382	K += kg_lt.K	×
383	energy_total += kg_lt.energy	×
384	energies["Contractility"] = kg_lt.energy	×
385
386	if Set.Contractility_external:	×
387	kg_ext = KgContractilityExternal(Geo)	×
388	kg_ext.compute_work(Geo, Set)	×
389	g += kg_ext.g	×
390	K += kg_ext.K	×
391	energy_total += kg_ext.energy	×
392	energies["Contractility_external"] = kg_ext.energy	×
393
394	# Substrate
395	if Set.Substrate == 2:	×
396	kg_subs = KgSubstrate(Geo)	×
397	kg_subs.compute_work(Geo, Set)	×
398	g += kg_subs.g	×
399	K += kg_subs.K	×
400	energy_total += kg_subs.energy	×
401	energies["Substrate"] = kg_subs.energy	×
402
403	# Triangle Energy Barrier
404	if Set.EnergyBarrierA:	×
405	kg_Tri = KgTriEnergyBarrier(Geo)	×
406	kg_Tri.compute_work(Geo, Set)	×
407	g += kg_Tri.g	×
408	K += kg_Tri.K	×
409	energy_total += kg_Tri.energy	×
410	energies["TriEnergyBarrier"] = kg_Tri.energy	×
411
412	# Triangle Energy Barrier Aspect Ratio
413	if Set.EnergyBarrierAR:	×
414	kg_TriAR = KgTriAREnergyBarrier(Geo)	×
415	kg_TriAR.compute_work(Geo, Set)	×
416	g += kg_TriAR.g	×
417	K += kg_TriAR.K	×
418	energy_total += kg_TriAR.energy	×
419	energies["TriEnergyBarrierAR"] = kg_TriAR.energy	×
420
421	if implicit_method is True:	×
422	# Viscous Energy
423	kg_Viscosity = KgViscosity(Geo)	×
424	kg_Viscosity.compute_work(Geo, Set, Geo_n)	×
425	g += kg_Viscosity.g	×
426	K += kg_Viscosity.K	×
427	energy_total = kg_Viscosity.energy	×
428	energies["Viscosity"] = kg_Viscosity.energy	×
429
430	return g, K, energy_total, energies	×
431
432
433	def gGlobal(Geo_0, Geo_n, Geo, Set, implicit_method=True):	1✔
434	"""
435	Compute the global gradient
436	:param Geo_0:
437	:param Geo_n:
438	:param Geo:
439	:param Set:
440	:param implicit_method:
441	:return:
442	"""
443	g = np.zeros((Geo.numY + Geo.numF + Geo.nCells) * 3, dtype=np.float64)	×
444	energies = {}	×
445
446	# Surface Energy
447	if Set.lambdaS1 > 0:	×
448	kg_SA = KgSurfaceCellBasedAdhesion(Geo)	×
449	kg_SA.compute_work(Geo, Set, None, False)	×
450	g += kg_SA.g	×
451	energies["Surface"] = kg_SA.energy	×
452
453	# Volume Energy
454	if Set.lambdaV > 0:	×
455	kg_Vol = KgVolume(Geo)	×
456	kg_Vol.compute_work(Geo, Set, None, False)	×
457	g += kg_Vol.g[:]	×
458	energies["Volume"] = kg_Vol.energy	×
459
460	if implicit_method is True:	×
461	# Viscous Energy
462	kg_Viscosity = KgViscosity(Geo)	×
463	kg_Viscosity.compute_work(Geo, Set, Geo_n, False)	×
464	g += kg_Viscosity.g	×
465	energies["Viscosity"] = kg_Viscosity.energy	×
466
467	# # TODO: Plane Elasticity
468	# if Set.InPlaneElasticity:
469	# gt, Kt, EBulk = KgBulk(geo_0, Geo, Set)
470	# K += Kt
471	# g += gt
472	# E += EBulk
473	# Energies["Bulk"] = EBulk
474
475	# Bending Energy
476	# TODO
477
478	# Triangle Energy Barrier
479	if Set.EnergyBarrierA:	×
480	kg_Tri = KgTriEnergyBarrier(Geo)	×
481	kg_Tri.compute_work(Geo, Set, None, False)	×
482	g += kg_Tri.g	×
483	energies["TriEnergyBarrier"] = kg_Tri.energy	×
484
485	# Triangle Energy Barrier Aspect Ratio
486	if Set.EnergyBarrierAR:	×
487	kg_TriAR = KgTriAREnergyBarrier(Geo)	×
488	kg_TriAR.compute_work(Geo, Set, None, False)	×
489	g += kg_TriAR.g	×
490	energies["TriEnergyBarrierAR"] = kg_TriAR.energy	×
491
492	# Propulsion Forces
493	# TODO
494
495	# Contractility
496	if Set.Contractility:	×
497	kg_lt = KgContractility(Geo)	×
498	kg_lt.compute_work(Geo, Set, None, False)	×
499	g += kg_lt.g	×
500	energies["Contractility"] = kg_lt.energy	×
501
502	# Contractility as external force
503	if Set.Contractility_external:	×
504	kg_ext = KgContractilityExternal(Geo)	×
505	kg_ext.compute_work(Geo, Set, None, False)	×
506	g += kg_ext.g	×
507	energies["Contractility_external"] = kg_ext.energy	×
508
509	# Substrate
510	if Set.Substrate == 2:	×
511	kg_subs = KgSubstrate(Geo)	×
512	kg_subs.compute_work(Geo, Set, None, False)	×
513	g += kg_subs.g	×
514	energies["Substrate"] = kg_subs.energy	×
515
516	return g, energies	×
517
518
519	def remeshing_cells(Geo_0, Geo_n, Geo, Dofs, Set, cells_to_change, ghost_node):	1✔
520	"""
521	Newton-Raphson method
522	:param Geo_0:
523	:param Geo_n:
524	:param Geo:
525	:param Dofs:
526	:param Set:s
527	:return:
528	"""
529
530	def objective_function(dy, dof):	×
531	# Reshape dy to match the geometry of the system
532	dy_reshaped = np.reshape(dy, (Geo.numF + Geo.numY + Geo.nCells, 3))	×
533
534	# Create a copy of Geo to not change the original one
535	Geo_copy = Geo.copy()	×
536
537	# Update the vertices in the copy
538	Geo_copy.update_vertices(dy_reshaped)	×
539
540	# Update the measures in the copy
541	Geo_copy.update_measures()	×
542
543	# Compute the energies
544	g, _, energy_total, _ = KgGlobal(Geo_0, Geo_n, Geo_copy, Set)	×
545
546	logger.info('Energy: ' + str(energy_total))	×
547	gr = np.linalg.norm(g[dof])	×
548	logger.info('\|\|gr\|\|= ' + str(gr))	×
549
550	# Return the total energy
551	return gr	×
552
553	def gradient_function(dy, dof):	×
554	# Reshape dy to match the geometry of the system
555	dy_reshaped = np.reshape(dy, (Geo.numF + Geo.numY + Geo.nCells, 3))	×
556
557	# Create a copy of Geo to not change the original one
558	Geo_copy = Geo.copy()	×
559
560	# Update the vertices in the copy
561	Geo_copy.update_vertices(dy_reshaped)	×
562
563	# Update the measures in the copy
564	Geo_copy.update_measures()	×
565
566	# Compute the energies
567	g, _, _, _ = KgGlobal(Geo_0, Geo_n, Geo_copy, Set)	×
568
569	g[dof == False] = 0	×
570
571	# Return the gradient
572	return g.flatten()	×
573
574	Geo.remeshing = True	×
575
576	if np.isin(ghost_node, Geo.XgBottom).any():	×
577	interface_type = 2	×
578	elif np.isin(ghost_node, Geo.XgTop).any():	×
579	interface_type = 0	×
580
581	Dofs.get_remeshing_dofs(Geo, cells_to_change, interface_type)	×
582
583	# Set the degrees of freedom
584	dof = Dofs.remeshing	×
585
586	dy_initial = np.zeros((Geo.numY + Geo.numF + Geo.nCells) * 3, dtype=np.float64)	×
587	g, K, _, _ = KgGlobal(Geo_0, Geo_n, Geo, Set)	×
588	dy_initial[dof] = ml_divide(K, dof, g)	×
589
590	lower_bounds = np.zeros_like(dy_initial)	×
591	upper_bounds = 10e-1 * np.ones_like(dy_initial)	×
592	bounds = Bounds(lower_bounds, upper_bounds)	×
593
594	result = minimize(	×
595	fun=objective_function, # Defined as before
596	x0=dy_initial,
597	args=(dof,),
598	method='L-BFGS-B', # Newton-CG
599	jac=gradient_function,
600	bounds=bounds,
601	options={'disp': True}
602	)
603
604	dy_reshaped = np.reshape(result.x, (Geo.numF + Geo.numY + Geo.nCells, 3))	×
605	return dy_reshaped	×

Pablo1990 / pyVertexModel / 11234788583

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