11836421490

Committed 14 Nov 2024 11:37AM UTC coverage: 0.789%. Remained the same

Build # 11836421490

Build Type

push

github

Committed by

Pablo1990

Commit Message

bugix

Run Details

0 of 1010 branches covered (0.0%)

Branch coverage included in aggregate %.

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

1 existing line in 1 file now uncovered.

51 of 5454 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

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
from src.pyVertexModel.geometry.face import get_interface

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:
    """
    # Bottom nodes
    all_bottom_nodes = []
    for cell in Geo.Cells:
        if cell.AliveStatus is not None:
            all_bottom_nodes.extend(cell.globalIds[np.any(np.isin(cell.T, Geo.XgBottom), axis=1)])
            for face in cell.Faces:
                if get_interface(face.InterfaceType) == get_interface('Bottom'):
                    all_bottom_nodes.append(face.globalIds)

    all_bottom_nodes_pos = all_bottom_nodes * 3
    all_bottom_nodes_pos.extend(3 * [i + 1 for i in all_bottom_nodes])
    all_bottom_nodes_pos.extend(3 * [i + 2 for i in all_bottom_nodes])
    dof_bottom = np.unique(all_bottom_nodes_pos)

    # Update the bottom nodes with the same displacement as the corresponding real nodes
    dy[dof, 0] = -Set.dt / Set.nu * g[dof]
    dy[dof_bottom, 0] = -Set.dt / Set.nu_bottom * g[dof_bottom]

    # 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]

                    # Check if the object has the attribute opposite_cell
                    if not hasattr(Geo.Cells[int(node)], 'opposite_cell'):
                        return dy

                    # 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: 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

Pablo1990 / pyVertexModel / 11836421490

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