utils_Quad4.py

import numpy as np
from numpy.polynomial.legendre import leggauss
import meshio

#-----------------------------------------------------------------------------#

def getPlaneStrain(Emod, nu):
    """
    Calculates the plane-strain stiffness matrix in Voigt notation.

    Args:
        Emod (float): Young's modulus.
        nu (float): Poisson's ratio.

    Returns:
       CMatx (np.ndarray): 4x4 stiffness matrix.
    """

    const = Emod*(1.0-nu)/((1+nu)*(1.0-2.0*nu))

    CMatx = np.zeros((3,3))
    CMatx[0,0] = const
    CMatx[1,1] = const
    CMatx[0,1] = const*nu/(1.0-nu)
    CMatx[1,0] = const*nu/(1.0-nu)
    CMatx[2,2] = (1.0-2.0*nu)*const/(2.0*(1.0-nu))
    
    return CMatx

#-----------------------------------------------------------------------------#

def gauss_2d(n):
    """
    Returns the integration point (natural) coordinates and weights 

    """
    
    xi_1d, wi_1d = leggauss(n)
    X, Y = np.meshgrid(xi_1d, xi_1d, indexing='ij')  
    W = np.outer(wi_1d, wi_1d)

    pts = np.column_stack([X.ravel(), Y.ravel()])   
    wts = W.ravel()                                 
    
    return pts, wts

#-----------------------------------------------------------------------------#

def getShapeFunc(xi, eta, nElNodes):
    """
    Calculates the shape functions for linear quad4 elements.

    Args:
        xi (float): xi natural coordinate of the integration point.
        eta (fload): eta natural coordinate of the integration point.
        nElNodes (int): number of nodes per element

    Returns:
        shapeFunc (ndarray): shape functions array
    """
    
    shapeFunc = np.zeros((1, nElNodes))

    shapeFunc[0][0] = (1.0-eta-xi+xi*eta)*0.25
    shapeFunc[0][1] = (1.0-eta+xi-xi*eta)*0.25
    shapeFunc[0][2] = (1.0+eta+xi+xi*eta)*0.25
    shapeFunc[0][3] = (1.0+eta-xi-xi*eta)*0.25
    
    return shapeFunc

#-----------------------------------------------------------------------------#

def getShapeFuncDeriv(xi, eta, nElNodes):
    """
    Calculate the shape function derivaties 

    Args:
        xi (float): xi natural coordinate of the integration point.
        eta (fload): eta natural coordinate of the integration point.
        nElNodes (int): number of nodes per element.
        
    Returns:
        shapeDeriv (ndarray): array of shape functions derivatives.
    """
    
    shapeDeriv = np.zeros((2, nElNodes))
        
    shapeDeriv[0,0] = (-1.0+eta)*0.25
    shapeDeriv[0,1] = (1.0-eta)*0.25
    shapeDeriv[0,2] = (1.0+eta)*0.25
    shapeDeriv[0,3] = (-1.0-eta)*0.25

    shapeDeriv[1,0] = (-1.0+xi)*0.25
    shapeDeriv[1,1] = (-1.0-xi)*0.25
    shapeDeriv[1,2] = (1.0+xi)*0.25
    shapeDeriv[1,3] = (1.0-xi)*0.25 
    
    return shapeDeriv

#-----------------------------------------------------------------------------#

def getElemDispDof(nElDispDofs, nElNodes, NodeConn):
    """
    Creates array of element displacement degrees of freedom.

    Args:
        nElDispDofs (int): number of element displacement dofs
        nElNodes (int): number of element nodes
        NodeConn (ndarray): array of nodal connectivity

    Returns:
        dispDof (ndarray): array of element displacement dofs
    """
    
    dispDof = np.zeros(nElDispDofs, dtype=int)
    
    for iNod in range(nElNodes):
        
        dispDof[2*iNod] = int(2*NodeConn[iNod])
        dispDof[2*iNod+1] = int(2*NodeConn[iNod]+1)
    
    return dispDof

#-----------------------------------------------------------------------------#

def getElemPhiDof(nElPhiDofs, nElNodes, NodeConn):
    """
    Creates array of element phi degrees of freedom.

    Args:
        nElPhiDofs (int): number of element phi dofs
        nElNodes (int): number of element nodes
        NodeConn (ndarray): array of nodal connectivity

    Returns:
        phiDof:  array of element phi dofs
    """
    
    phiDof = np.zeros(nElPhiDofs, dtype=int)
    
    for iNod in range(nElNodes):
        
        phiDof[iNod] = int(NodeConn[iNod])
    
    return phiDof

#-----------------------------------------------------------------------------#

def CalcCartDeriv(elNodCoord, nElNodes, nStres, nElDispDofs, sFuncDeriv, wt):
    """
    Calculates the derivative matrix BMAT, strain matrix BuMat and integration point volume intVol.

    Args:
        elNodCoord (_type_): _description_
        nElNodes (_type_): _description_
        nStres (_type_): _description_
        nElDispDofs (_type_): _description_
        sFuncDeriv (_type_): _description_
        wt (_type_): _description_

    Returns:
        _type_: _description_
    """
    
    jacMat = np.dot(sFuncDeriv, elNodCoord)
    
    intVol = np.linalg.det(jacMat)*wt
    
    cartDeriv = np.dot(np.linalg.inv(jacMat), sFuncDeriv)
    
    strainMat = np.zeros((nStres, nElDispDofs))
    
    for iNod in range(nElNodes):
        
        strainMat[0,2*iNod]   = cartDeriv[0,iNod];
        strainMat[0,2*iNod+1] = 0;
        strainMat[1,2*iNod]   = 0;
        strainMat[1,2*iNod+1] = cartDeriv[1,iNod];
        strainMat[2,2*iNod]   = cartDeriv[1,iNod];
        strainMat[2,2*iNod+1] = cartDeriv[0,iNod];
        
    return intVol, cartDeriv, strainMat

#-----------------------------------------------------------------------------#

def CalcElStran(nElements, nGauss, nStres, elemDispDof, BuMat, u):
    """
    Calculates the strain tensro in voigt notation at integratin points from the displacement vector.

    Args:
        nElements (_type_): _description_
        nGauss (_type_): _description_
        nStres (_type_): _description_
        elemDispDof (_type_): _description_
        BuMat (_type_): _description_
        u (_type_): _description_

    Returns:
        _type_: _description_
    """
    
    elStran = np.zeros((nElements, nGauss, nStres))
   
    for iElem in range(nElements):
        
        elemDOFs = elemDispDof[iElem]
        
        for iGaus in range(nGauss):
            
            xx =  u[elemDOFs]
            
            elStran[iElem,iGaus] = np.dot(BuMat[iElem,iGaus], u[elemDOFs])
            
    return elStran

#-----------------------------------------------------------------------------#

def Mises(S, nu):
    """
    Von-Mises equivalent stress in plane strain
    """
    sx = S[0]
    sy = S[1]
    txy = S[2]
    sz = nu*(sx+sy)
    
    return np.sqrt( 0.5 * ( (sx - sy)**2 + (sy - sz)**2 + (sz-sx)**2 + 3*txy**2 ) )

#-----------------------------------------------------------------------------#

def R_pow(k, n, p):    # Power-law hardening    
    return k*p**n

#-----------------------------------------------------------------------------#

def dR_pow(k, n, p):    
    return k*n*p**(n-1)

#-----------------------------------------------------------------------------#

def CalcFint(elStres, nElements, nGauss, elemDispDof, nTotDof, BuMat, intPtVol):
    """
    Calculates the internal force vector

    Args:
        elStres (_type_): _description_
        nElements (_type_): _description_
        nGauss (_type_): _description_
        elemDispDof (_type_): _description_
        nTotDof (_type_): _description_
        BuMat (_type_): _description_
        intPtVol (_type_): _description_

    Returns:
        _type_: _description_
    """
    
    Fint = np.zeros((nTotDof))

    for iElem in range(nElements):

        elemDOFs = elemDispDof[iElem]
        
        for iGaus in range(nGauss):
                       
            Fint[elemDOFs] += np.dot(BuMat[iElem,iGaus].T, elStres[iElem,iGaus])*intPtVol[iElem,iGaus]
        
    return Fint

#-----------------------------------------------------------------------------#

def CalcNodValsPlastic(nElements, nGauss, nStres, nNodes, elemNodeConn, elStres, elStran, elStran_e, elStran_p, elStres_eq, elStran_eq, intPtVol):
    """
    Maps the integratin point values like stress and strain to the nodes. This is required for visualization in ParaView.

    Args:
        nElements (_type_): _description_
        nGauss (_type_): _description_
        nStres (_type_): _description_
        nNodes (_type_): _description_
        elemNodeConn (_type_): _description_
        elStres (_type_): _description_
        elStran (_type_): _description_
        elStran_e (_type_): _description_
        elStran_p (_type_): _description_
        elStres_eq (_type_): _description_
        elStran_eq (_type_): _description_
        intPtVol (_type_): _description_

    Returns:
        _type_: _description_
    """
    
    nodStran = np.zeros((nNodes, nStres))
    nodStran_e = np.zeros((nNodes, nStres))
    nodStran_p = np.zeros((nNodes, nStres))
    nodStres = np.zeros((nNodes, nStres))
    nodStran_eq = np.zeros((nNodes))
    nodStres_eq = np.zeros((nNodes))
    nodCountVol = np.zeros(nNodes)

    for iElem in range(nElements):
                
        for iGaus in range(nGauss):
            
            for iNod in elemNodeConn[iElem]:
                
                nodStran[iNod] += elStran[iElem,iGaus]*intPtVol[iElem,iGaus]
                nodStran_e[iNod] += elStran_e[iElem,iGaus]*intPtVol[iElem,iGaus]
                nodStran_p[iNod] += elStran_p[iElem,iGaus]*intPtVol[iElem,iGaus]
                nodStres[iNod] += elStres[iElem,iGaus]*intPtVol[iElem,iGaus]
                nodStran_eq[iNod] += elStran_eq[iElem,iGaus]*intPtVol[iElem,iGaus]
                nodStres_eq[iNod] += elStres_eq[iElem,iGaus]*intPtVol[iElem,iGaus]
                nodCountVol[iNod] += intPtVol[iElem,iGaus]
                
    for iNod in range(nNodes):
    
        nodStres[iNod] = nodStres[iNod]/nodCountVol[iNod]
        nodStran[iNod] = nodStran[iNod]/nodCountVol[iNod]
        nodStran_e[iNod] = nodStran_e[iNod]/nodCountVol[iNod]
        nodStran_p[iNod] = nodStran_p[iNod]/nodCountVol[iNod]
        nodStran_eq[iNod] = nodStran_eq[iNod]/nodCountVol[iNod]
        nodStres_eq[iNod] = nodStres_eq[iNod]/nodCountVol[iNod]
        
    return nodStres, nodStran, nodStran_e, nodStran_p, nodStres_eq, nodStran_eq 

#-----------------------------------------------------------------------------#

def RMIsoHard_PFF(EG3, ENU, CMatx_e, K_hard, n_pow, sig_y0, iStep, elDStran, elStran_e, elStran_p, elStres, elStran_eq, I, gPhi_d, w_pOld):
    """
    Return mapping algorithm modified for phase-field fracture.

    Args:
        EG3 (_type_): _description_
        ENU (_type_): _description_
        CMatx_e (_type_): _description_
        K_hard (_type_): _description_
        n_pow (_type_): _description_
        sig_y0 (_type_): _description_
        iStep (_type_): _description_
        elDStran (_type_): _description_
        elStran_e (_type_): _description_
        elStran_p (_type_): _description_
        elStres (_type_): _description_
        elStran_eq (_type_): _description_
        I (_type_): _description_
        gPhi_d (_type_): _description_
        w_pOld (_type_): _description_

    Returns:
        _type_: _description_

    """
    
    # Elastic strain
    elStran_e[:] += elDStran
    
    # Trial stress
    sig_trial = np.dot(CMatx_e, elStran_e) * gPhi_d
    
    # Von Mises stress
    sig_trial_eq = Mises(sig_trial, ENU)
    
    # Yield condition
    deqpl = 2.220446049250313e-16     # Plastic strain increment
    p = elStran_eq + deqpl            # Total equivalent plastic strain
    
    sYield0 = sig_y0 + R_pow(K_hard, n_pow, p) # Current yield stress
    hard = dR_pow(K_hard, n_pow, p)            # Hardening modulus
    
    f_yield = sig_trial_eq - sYield0*(1.0 + 1.0e-6) # Yield function
                    
    # Check yielding
    if f_yield <= 0:    # --> Elastic step
        
        elStres[:] = sig_trial
        
        return sig_trial_eq, elStran_eq, w_pOld, CMatx_e*gPhi_d

    else:    # --> Plastic step
        
        sig_trial_dev = sig_trial - 1/3*(sig_trial[0]+sig_trial[1]+ENU*(sig_trial[0]+sig_trial[1]))*I  # Deviatoric stress
        N_tr = 3/2 * sig_trial_dev/sig_trial_eq   # Plastic flow direction

        nIter_RM = 0    # Iteration counter
    
        # Return mapping algorithm
        while abs(f_yield) > 1e-3:
            
            # Update Iteration counter
            nIter_RM += 1
            
            # Update plastic strain increment --> Power-law hardening
            deqpl += f_yield/(EG3 * gPhi_d + hard)  # Plastic strain increment
            p = elStran_eq + deqpl
            sYield = sig_y0 + R_pow(K_hard, n_pow, p)
            hard = dR_pow(K_hard, n_pow, p)
            f_yield = sig_trial_eq - EG3 * gPhi_d * deqpl - sYield # Update yield condition
                        
            if nIter_RM == 20:
                print(f'Return mapping did not converge at step {iStep} 😭😭😭')
                break
        
        # ----- CONVERGED 🥳🥳🥳 -----
        
        dep = deqpl*N_tr
        elStran_p[:] += dep                   # Plastic strain
        elStran_e[:] -= dep                   # Elastic strain 
        
        elStres[:] = np.dot(CMatx_e, elStran_e)      # Stress
                
        w_p = elStres[0]*dep[0] + elStres[1]*dep[1] + 2*elStres[2]*dep[2]
        w_p += w_pOld
                
        elStres[:] = elStres[:] * gPhi_d # Stress
        
        sig_eq = Mises(elStres, ENU)     # Equivalent stress
        eps_eq = p                       # Equivalent plastic strain
                
        Ce_N = np.dot(CMatx_e, N_tr)
        N_Ce_N = np.dot(Ce_N, N_tr)
        CMatx = (CMatx_e - (np.outer(Ce_N, Ce_N) / (2/3*hard + N_Ce_N) )) * gPhi_d
    
        return sig_eq, eps_eq, w_p, CMatx

#-----------------------------------------------------------------------------#

def voigt2tensor(E_voigt):
    """
    Converts strain from voigt to tensorial notation required for spectral decomposition. 

    Args:
        E_voigt (_type_): _description_

    Returns:
        _type_: _description_
    """
    E_tensor = np.zeros((2, 2))
    E_tensor[0, 0] = E_voigt[0]
    E_tensor[1, 1] = E_voigt[1]
    E_tensor[0, 1] = E_tensor[1, 0] = E_voigt[2] / 2
    return E_tensor

#-----------------------------------------------------------------------------#

def CalcPsiSpectral(nElements, nGauss, lam, G, elStran_e, psi_plus, psi_minus):
    """
    Performs spectral split to the elastic strain tensor and calculates the 
    positive and negative parts of the elsastic strain energy density according to 
    Miehe et al. (2010).

    Args:
        nElements (_type_): _description_
        nGauss (_type_): _description_
        lam (_type_): _description_
        G (_type_): _description_
        elStran_e (_type_): _description_
        psi_plus (_type_): _description_
        psi_minus (_type_): _description_
    """
    
    for iElem in range(nElements):
                
        for iGauss in range(nGauss):
                        
            # Strain tensor
            stran_tensor = voigt2tensor(elStran_e[iElem, iGauss])
            
            # Volumetric parts (with Macauley brackets)
            tr_eps = np.trace(stran_tensor)
            tr_plus = max(tr_eps, 0.0)
            tr_minus = min(tr_eps, 0.0)
            
            # Eigen values
            eigvals, eigvecs = np.linalg.eigh(stran_tensor)
            
            # Positive/Negative parts of the strain
            E_plus = sum(np.maximum(e, 0) * np.outer(v, v) for e, v in zip(eigvals, eigvecs.T))
            E_minus = sum(np.minimum(e, 0) * np.outer(v, v) for e, v in zip(eigvals, eigvecs.T))
            
            # Positive/Negative parts of the strain energy density
            psi_plus[iElem, iGauss] = 0.5*lam*tr_plus**2 +  G*np.trace(np.dot(E_plus , E_plus))
            psi_minus[iElem, iGauss] = 0.5*lam*tr_minus**2 +  G*np.trace(np.dot(E_minus , E_minus))

#-----------------------------------------------------------------------------#

def CalcDrivForcBrittle(nElements, nGauss, psi_plus, wc, elemH):
    """
    Calculates the crack brittle driving force. 

    Args:
        nElements (_type_): _description_
        nGauss (_type_): _description_
        psi_plus (_type_): _description_
        wc (_type_): _description_
        elemH (_type_): _description_
    """
    
    for iElem in range(nElements):
        for iGauss in range(nGauss):
            
            drivForce = (psi_plus[iElem, iGauss])/wc
            elemH[iElem, iGauss] = max(drivForce, elemH[iElem, iGauss])
            
#-----------------------------------------------------------------------------#

def CalcDrivForcEP(nElements, nGauss, psi_plus, wc, w_p, elemH):
    """
    Calculates the elastoplastic crack driving force. 

    Args:
        nElements (_type_): _description_
        nGauss (_type_): _description_
        psi_plus (_type_): _description_
        wc (_type_): _description_
        w_p (_type_): _description_
        elemH (_type_): _description_
    """
    
    for iElem in range(nElements):
        for iGauss in range(nGauss):
            
            drivForce = (psi_plus[iElem, iGauss] + w_p[iElem, iGauss])/wc
            elemH[iElem, iGauss] = max(drivForce, elemH[iElem, iGauss])

#-----------------------------------------------------------------------------#

def CalcDrivForcEP_TH(nElements, nGauss, psi_plus, wc, w_p, elemH):
    """
    Calculates the elastoplastic crack driving force with a threshold. 

    Args:
        nElements (_type_): _description_
        nGauss (_type_): _description_
        psi_plus (_type_): _description_
        wc (_type_): _description_
        w_p (_type_): _description_
        elemH (_type_): _description_
    """
    
    for iElem in range(nElements):
        for iGauss in range(nGauss):
            
            drivForce = max((psi_plus[iElem, iGauss] + w_p[iElem, iGauss])/wc - 1, 0)
            elemH[iElem, iGauss] = max(drivForce, elemH[iElem, iGauss])

#-----------------------------------------------------------------------------#

def CalcElPhi(nElements, nGauss, elemPhiDof, shapeFunc, phi):
    """
    Calculates phi values at integration points. 

    Args:
        nElements (_type_): _description_
        nGauss (_type_): _description_
        elemPhiDof (_type_): _description_
        shapeFunc (_type_): _description_
        phi (_type_): _description_

    Returns:
        _type_: _description_
    """
    
    elPhi = np.zeros((nElements, nGauss))
    
    for iElem in range(nElements):
        
        for iGauss in range(nGauss):
            
            elPhi[iElem, iGauss] = np.dot(shapeFunc[iGauss], phi[elemPhiDof[iElem]])[0]
            
    return elPhi

#-----------------------------------------------------------------------------#

def ConstrainElPhi(nElements, nGauss, elPhi, psi_plus, psi_minus):
    """
    Constrains phi according to Eq.(27c) in Ambati et al(2015)

    Args:
        nElements (float): Number of elements.
        nGauss (float): Number of integration points.
        elPhi (ndarray): Integration point phi. 
        psi_plus (ndarray): Positive part of the strain energy.
        psi_minus (ndarray): Negative part of the strain energy.
    """
        
    for iElem in range(nElements):
        
        for iGauss in range(nGauss):
            
            if psi_minus[iElem, iGauss] > psi_plus[iElem, iGauss]:
                elPhi[iElem, iGauss] = 0

#-----------------------------------------------------------------------------#

def CalcFp(nElements, nGauss, nTotPhiDof, elemPhiDof, shapeFunc, intPtVol, elemH):
    """
    Calculates the right hand side of the phase-field fracture system of equation. 

    Args:
        nElements (_type_): _description_
        nGauss (_type_): _description_
        nTotPhiDof (_type_): _description_
        elemPhiDof (_type_): _description_
        shapeFunc (_type_): _description_
        intPtVol (_type_): _description_
        elemH (_type_): _description_

    Returns:
        _type_: _description_
    """
    
    Fp = np.zeros((nTotPhiDof))
    
    for iElem in range(nElements):
        
        elemPhiDOFs = elemPhiDof[iElem]
        
        for iGauss in range(nGauss):
            
            dummyVar = shapeFunc[iGauss].T*elemH[iElem, iGauss]*intPtVol[iElem, iGauss]
                                    
            Fp[elemPhiDOFs] += dummyVar.flatten() 
        
    return Fp

#-----------------------------------------------------------------------------#

def getTopY_DOFs(mesh, ly, nDim):
    """
    Gets the y dofs for top nodes. 

    Args:
        mesh (_type_): _description_
        ly (_type_): _description_
        nDim (_type_): _description_

    Returns:
        _type_: _description_
    """
    
    topY_DOFs = []

    # Number of nodes
    nNodes = mesh.points.shape[0]  
    # Node coordinates
    nodeCoord = mesh.points

    for iNod in range(nNodes):
        # top nodes
        if nodeCoord[iNod][1]==ly:
            topY_DOFs.append(nDim*iNod + 1)
    
    return topY_DOFs

#-----------------------------------------------------------------------------#

def TensileDisp2D(ly, yDisp, mesh):
    """
    Applies tensile displacement boundary conditions to a regular quadrilateral in the y direction. The origin point must be (0,0).

    Args:
        ly (float): y-length
        yDisp (float): y-displacement
        mesh (meshio): Mesh object
    Returns:
        list: List containing [[nod id, dof, value]]
    """
    
    # List of prescribed degrees of freedom. Order of list [node id, dof, value]
    presBCs = []   
    # Number of nodes
    nNodes = mesh.points.shape[0]  
    # Node coordinates
    nodeCoord = mesh.points[:,0:2]
    
    # Loop through nodes  
    for iNod in range(nNodes):
        # Bottom nodes
        if nodeCoord[iNod][1]==0:
            # Find bottom left corner node.
            if nodeCoord[iNod][0] == 0: 
                # Apply fixed BC
                presBCs.append([iNod, 0, 0])
                presBCs.append([iNod, 1, 0]) 
            # Other bottom nodes
            else :  
                # Y-fixed   
                presBCs.append([iNod, 1, 0])

        # Top nodes
        elif nodeCoord[iNod][1]==ly:
            # Top left corner node
            if nodeCoord[iNod][0] == 0:
                # Fix-x 
                presBCs.append([iNod, 0, 0])
                presBCs.append([iNod, 1, yDisp])
            # Other top nodes
            else :
                # Prescribed y
                presBCs.append([iNod, 1, yDisp])
    
    return presBCs
    
#-----------------------------------------------------------------------------#

def WriteDispBCs(Simul, elementName, mesh, presBCs, dispDofs=2):
    """
    Function to help visualize the displacement BCs in Paraview. Writes mesh with prescibed BCs as vtk.

    Args:
        Simul (str): Simulation name
        elementName (str): meshio compatible element name.
        mesh (meshio): Mesh object.
        presBCs (list): List containing [[nod id, dof, value]].
        dispDofs (int, optional): Displacement dofs. Defaults to 2.

    Returns:
        meshio: Mesh object with BCs.
    """
       
    # Number of prescribed dofs
    nPresDofs = len(presBCs)
    # Number of nodes
    nNodes = mesh.points.shape[0] 
    # Node coordinates
    nodeCoord = mesh.points
    # Node connectivity
    nodeConnectivity = mesh.cells_dict[elementName] 
    
    # Vector of displacement dofs values (like solution vector)
    uDisp = np.zeros(dispDofs*nNodes)
    # Vector of displacement dofs flag 
    flagsDisp = np.zeros(dispDofs*nNodes, dtype=int)
    
    # Vector of displacement dofs values (compatible with meshio)
    sdisp = np.zeros((nNodes, dispDofs))
    # Vector of displacement dofs flag
    fdisp = np.zeros((nNodes, dispDofs), dtype=int)

    # Create new mesh object for writing BCs  
    cells = [
        (elementName, nodeConnectivity),
    ]
    BCmesh = meshio.Mesh(
        nodeCoord,
        cells,
    )
    
    # Arrange BCs as solution vector  
    for i in range(nPresDofs):
        pDOF = presBCs[i][0]*dispDofs + presBCs[i][1]
        uDisp[pDOF] = presBCs[i][2]
        flagsDisp[pDOF] = 1
        
    # Rearrange BCs array of `nNodes x dispDofs` for meshio output
    n = nNodes*dispDofs

    sdisp[:,0] = uDisp[0:n:dispDofs]
    sdisp[:,1] = uDisp[1:n:dispDofs]         

    fdisp[:,0] = flagsDisp[0:n:dispDofs]
    fdisp[:,1] = flagsDisp[1:n:dispDofs]
    
    if dispDofs == 3:
        sdisp[:,2] = uDisp[2:n:dispDofs]
        fdisp[:,2] = flagsDisp[2:n:dispDofs]
    
    
    # Append and write to vtk
    BCmesh.point_data.update({"disp": sdisp})
    BCmesh.point_data.update({"FlagBC": fdisp})
    BCmesh.write(Simul+"_BC.vtu")
    
    return BCmesh

#-----------------------------------------------------------------------------#