11251700108

Committed 09 Oct 2024 09:01AM UTC coverage: 8.796% (+0.007%) from 8.789%

Build # 11251700108

Build Type

push

github

Committed by

Pablo1990

Commit Message

In case there are no periodic boundaries

Run Details

0 of 1004 branches covered (0.0%)

Branch coverage included in aggregate %.

1 of 22 new or added lines in 1 file covered. (4.55%)

2 existing lines in 1 file now uncovered.

564 of 5408 relevant lines covered (10.43%)

0.1 hits per line

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

6.25

/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
    dy = map_vertices_periodic_boundaries(Geo, dy)

    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 map_vertices_periodic_boundaries(Geo, dy):
    """
    Update the border ghost nodes with the same displacement as the corresponding real nodes
    :param Geo:
    :param dy:
    :return:
    """
    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
                    if Geo.Cells[int(numCell)].opposite_cell is None:
                        return dy
                    
                    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: what happens if it is an scutoid
                        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]

                            # NOW IT IS 0, WITH THE AVERAGE IT THERE WERE ERRORS
                            dy[global_id * 3:global_id * 3 + 3, 0] = 0
                        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]
    return dy


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

Pablo1990 / pyVertexModel / 11251700108

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