Skip to content

Commit

Permalink
Squashed commit of the following:
Browse files Browse the repository at this point in the history 
commit 83db852
Author: matnoel <[email protected]>
Date:   Wed Aug 14 09:48:20 2024 +0200

    Update pyproject.toml

commit ede6a75
Author: matnoel <[email protected]>
Date:   Wed Aug 14 09:48:00 2024 +0200

    Update CHANGELOG.md

commit 146cc8b
Author: matnoel <[email protected]>
Date:   Wed Aug 14 09:46:05 2024 +0200

    Update PostProcess.py

commit 2e3d656
Author: matnoel <[email protected]>
Date:   Wed Aug 14 09:45:33 2024 +0200

    Update phasefield examples

commit 0b11aeb
Author: matnoel <[email protected]>
Date:   Wed Aug 14 09:38:50 2024 +0200

    Update _dic.py

commit 50760a1
Author: matnoel <[email protected]>
Date:   Wed Aug 14 09:31:27 2024 +0200

    Update _phasefield.py

commit 8a7728a
Author: matnoel <[email protected]>
Date:   Wed Aug 14 09:20:01 2024 +0200

    Update _dic.py

commit bdb2b95
Author: matnoel <[email protected]>
Date:   Tue Aug 13 18:44:01 2024 +0200

    Update _dic.py

commit 04547d1
Author: matnoel <[email protected]>
Date:   Mon Aug 12 18:56:44 2024 +0200

    Update _gmsh_interface.py

commit 826b1e5
Author: matnoel <[email protected]>
Date:   Mon Aug 12 11:13:56 2024 +0200

    Update _dic.py

commit b7dac81
Author: matnoel <[email protected]>
Date:   Mon Aug 12 09:07:33 2024 +0200

    Update _elastic.py

commit 96f9f77
Author: matnoel <[email protected]>
Date:   Mon Aug 12 09:07:24 2024 +0200

    Update Geoms.py

commit 76f4e9d
Author: matnoel <[email protected]>
Date:   Wed Aug 7 21:03:17 2024 +0200

    Update _dic.py

    Clarifies implementation to be consistent with thesis manuscript.

commit 9f71c9c
Author: matnoel <[email protected]>
Date:   Tue Aug 6 17:41:55 2024 +0200

    Update _dic.py

commit 30af24d
Author: matnoel <[email protected]>
Date:   Mon Aug 5 19:17:45 2024 +0200

    Update _dic.py

commit b9861fc
Author: matnoel <[email protected]>
Date:   Mon Aug 5 18:33:22 2024 +0200

    Update _group_elems.py

commit 6e8cdc0
Author: matnoel <[email protected]>
Date:   Mon Aug 5 17:40:25 2024 +0200

    Update _group_elems.py

commit 934e81d
Author: matnoel <[email protected]>
Date:   Mon Aug 5 17:34:49 2024 +0200

    Update _mesh.py

commit 62a78b9
Author: matnoel <[email protected]>
Date:   Mon Aug 5 16:36:14 2024 +0200

    Update _group_elems.py

commit e524932
Author: matnoel <[email protected]>
Date:   Mon Aug 5 16:35:30 2024 +0200

    Update PyVista_Interface.py

commit 1780444
Author: matnoel <[email protected]>
Date:   Tue Jul 30 19:18:21 2024 +0200

    Update _dic.py

commit f999a83
Author: matnoel <[email protected]>
Date:   Tue Jul 30 19:17:20 2024 +0200

    Update _elastic.py

commit d6e29a4
Author: matnoel <[email protected]>
Date:   Tue Jul 30 09:33:04 2024 +0200

    Update _phasefield.py

commit c02b910
Author: matnoel <[email protected]>
Date:   Tue Jul 30 09:32:23 2024 +0200

    Update Display.py

commit c41cd57
Author: matnoel <[email protected]>
Date:   Tue Jul 30 09:32:11 2024 +0200

    Update Geoms.py

commit e4c5d1d
Author: matnoel <[email protected]>
Date:   Tue Jul 16 19:13:57 2024 +0200

    Update Solvers.py

commit 4d7f260
Author: matnoel <[email protected]>
Date:   Tue Jul 16 09:44:58 2024 +0200

    Update PhaseField

commit d2487af
Author: matnoel <[email protected]>
Date:   Wed Jul 10 15:06:07 2024 +0200

    Update Materials_test.py

commit 883b3fa
Merge: cbb79da ec5d437
Author: matnoel <[email protected]>
Date:   Wed Jul 3 12:59:48 2024 +0200

    Merge branch 'dev' of https://github.com/matnoel/EasyFEA into dev

commit cbb79da
Author: matnoel <[email protected]>
Date:   Tue Jul 2 10:25:27 2024 +0200

    Update simu.Set_Iter()

    In the previous version, activating Set_Iter() triggered the mandatory updating of the history field in damage simulations. This resulted in post-processing times that were far too long for the creation of videos or paraviews.

commit ec5d437
Author: matnoel <[email protected]>
Date:   Tue Jul 2 10:25:27 2024 +0200

    Update simu.Set_Iter()

    In the previous version, activating Set_Iter() triggered the mandatory updating of the history field in damage simulations. This resulted in post-processing times that were far too long for the creation of videos or paraviews.
  • Loading branch information
@matnoel
matnoel committed Aug 14, 2024
1 parent 78e34da commit 886852f
Show file tree
Hide file tree
Showing 22 changed files with 639 additions and 423 deletions.
11 changes: 11 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,17 @@

This document describes the changes made to the project.

## 1.2.0 (August 14, 2024):

- Added the resetAll option in Set_Iter() to simplify the update process after iteration activation.
- Enhanced clarity in phase field functions, including both simulation and material aspects.
- Improved display options for geometric objects.
- Improved display functions.
- Provided clearer functionality in mesh and group element.
- Updated the interface with Gmsh.
- Improved phasefield examples.
- Enhanced the Digital Image Correlation (DIC) analysis module.

## 1.1.0 (June 29, 2024):

- Updated gmsh interface functions.
Expand Down
6 changes: 3 additions & 3 deletions EasyFEA/Geoms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -293,7 +293,7 @@ def Symmetry(self, point=(0,0,0), n=(1,0,0)) -> None:
dec = newCoord - oldCoord
[point.Translate(*dec[p]) for p, point in enumerate(self.points)]

def Plot(self, ax: plt.Axes=None, color:str="", name:str="", plotPoints=True) -> plt.Axes:
def Plot(self, ax: plt.Axes=None, color:str="", name:str="", lw=None, ls=None, plotPoints=True) -> plt.Axes:

from EasyFEA import Display

Expand All @@ -304,9 +304,9 @@ def Plot(self, ax: plt.Axes=None, color:str="", name:str="", plotPoints=True) ->

lines, points = self.Get_coord_for_plot()
if color != "":
ax.plot(*lines[:,:inDim].T, color=color, label=self.name)
ax.plot(*lines[:,:inDim].T, color=color, label=self.name, lw=lw, ls=ls)
else:
ax.plot(*lines[:,:inDim].T, label=self.name)
ax.plot(*lines[:,:inDim].T, label=self.name, lw=lw, ls=ls)
if plotPoints:
ax.plot(*points[:,:inDim].T, ls='', marker='.',c='black')

Expand Down
15 changes: 11 additions & 4 deletions EasyFEA/fem/_gmsh_interface.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1035,10 +1035,14 @@ def Mesh_2D(self, contour: _Geom, inclusions: list[_Geom]=[], elemType=ElemType.
# Crack creation
crackLines, crackSurfaces, openPoints, openLines = self._Cracks_SetPhysicalGroups(cracks, entities2D)

if (len(cracks) > 0 and 'QUAD' in elemType) or len(additionalSurfaces) > 0:
# dont delete
surfaceTags = [s[1] for s in gmsh.model.getEntities(2)]
self._OrganiseSurfaces(surfaceTags, elemType, isOrganised)
# get created surfaces
surfaces = [entity[1] for entity in factory.getEntities(2)]
self._OrganiseSurfaces(surfaces, elemType, isOrganised)
# if (len(cracks) > 0 and 'QUAD' in elemType) or len(additionalSurfaces) > 0:
# # dont delete
# self._Synchronize()
# surfaceTags = [s[1] for s in gmsh.model.getEntities(2)]
# self._OrganiseSurfaces(surfaceTags, elemType, isOrganised)

self._Refine_Mesh(refineGeoms, contour.meshSize)

Expand Down Expand Up @@ -1721,6 +1725,9 @@ def types(elemType: str):

dict_results: dict[str, list[np.ndarray]] = {result: [] for result in results}

# activates the first iteration
simu.Set_Iter(0, resetAll=True)

for i in range(simu.Niter):
simu.Set_Iter(i)
[dict_results[result].append(reshape(simu.Result(result))) for result in results]
Expand Down
106 changes: 41 additions & 65 deletions EasyFEA/fem/_group_elems.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ def coordGlob(self) -> np.ndarray:
@coordGlob.setter
def coordGlob(self, coord: np.ndarray) -> None:
if coord.shape == self.__coordGlob.shape:
self.__coordGlob = coord
self.__coordGlob = coord
self._InitMatrix()

@property
Expand Down Expand Up @@ -308,7 +308,7 @@ def Get_dN_e_pg(self, matrixType: MatrixType) -> np.ndarray:
dN_pg = self.Get_dN_pg(matrixType)

# Derivation of shape functions in the real base
dN_e_pg: np.ndarray = np.einsum('epik,pkj->epij', invF_e_pg, dN_pg, optimize='optimal')
dN_e_pg: np.ndarray = np.einsum('epdk,pkn->epdn', invF_e_pg, dN_pg, optimize='optimal')
self.__dict_dN_e_pg[matrixType] = dN_e_pg

return self.__dict_dN_e_pg[matrixType].copy()
Expand All @@ -329,7 +329,7 @@ def Get_ddN_e_pg(self, matrixType: MatrixType) -> np.ndarray:

ddN_pg = self.Get_ddN_pg(matrixType)

ddN_e_pg = np.array(np.einsum('epik,pkj->epij', invF_e_pg, ddN_pg, optimize='optimal'))
ddN_e_pg = np.array(np.einsum('epdk,pkn->epdn', invF_e_pg, ddN_pg, optimize='optimal'))
self.__dict_ddN_e_pg[matrixType] = ddN_e_pg

return self.__dict_ddN_e_pg[matrixType].copy()
Expand All @@ -345,11 +345,7 @@ def Get_Nv_e_pg(self) -> np.ndarray:

invF_e_pg = self.Get_invF_e_pg(matrixType)
Nv_pg = self.Get_Nv_pg(matrixType)

Ne = self.Ne
nPe = self.nPe
jacobian_e_pg = self.Get_jacobian_e_pg(matrixType)
pg = self.Get_gauss(matrixType)

Nv_e_pg = invF_e_pg @ Nv_pg

Expand All @@ -372,11 +368,7 @@ def Get_dNv_e_pg(self) -> np.ndarray:

invF_e_pg = self.Get_invF_e_pg(matrixType)
dNv_pg = self.Get_dNv_pg(matrixType)

Ne = self.Ne
nPe = self.nPe
jacobian_e_pg = self.Get_jacobian_e_pg(matrixType)
pg = self.Get_gauss(matrixType)

dNv_e_pg = invF_e_pg @ dNv_pg

Expand All @@ -398,12 +390,8 @@ def Get_ddNv_e_pg(self) -> np.ndarray:
matrixType = MatrixType.beam

invF_e_pg = self.Get_invF_e_pg(matrixType)
ddNv_pg = self.Get_ddNv_pg(matrixType)

Ne = self.Ne
ddNv_pg = self.Get_ddNv_pg(matrixType)
nPe = self.nPe
jacobian_e_pg = self.Get_jacobian_e_pg(matrixType)
pg = self.Get_gauss(matrixType)

ddNv_e_pg = invF_e_pg @ invF_e_pg @ ddNv_pg

Expand Down Expand Up @@ -490,7 +478,7 @@ def Get_leftDispPart(self, matrixType: MatrixType) -> np.ndarray:

def Get_ReactionPart_e_pg(self, matrixType: MatrixType) -> np.ndarray:
"""Returns the part that builds the reaction term (scalar).
ReactionPart_e_pg = jacobian_e_pg * weight_pg * r_e_pg * N_pg' * N_pg\n
ReactionPart_e_pg = r_e_pg * jacobian_e_pg * weight_pg * N_pg' * N_pg\n
Returns -> jacobian_e_pg * weight_pg * N_pg' * N_pg
"""
Expand All @@ -509,32 +497,30 @@ def Get_ReactionPart_e_pg(self, matrixType: MatrixType) -> np.ndarray:

return self.__dict_phaseField_ReactionPart_e_pg[matrixType].copy()

def Get_DiffusePart_e_pg(self, matrixType: MatrixType, A: np.ndarray) -> np.ndarray:
def Get_DiffusePart_e_pg(self, matrixType: MatrixType) -> np.ndarray:
"""Returns the part that builds the diffusion term (scalar).
DiffusePart_e_pg = jacobian_e_pg * weight_pg * k * dN_e_pg' * A * dN_e_pg\n
DiffusePart_e_pg = k_e_pg * jacobian_e_pg * weight_pg * dN_e_pg' * A * dN_e_pg\n
Returns -> jacobian_e_pg * weight_pg * dN_e_pg' * A * dN_e_pg
Returns -> jacobian_e_pg * weight_pg * dN_e_pg'
"""

assert matrixType in MatrixType.Get_types()

if matrixType not in self.__dict_DiffusePart_e_pg.keys():

assert len(A.shape) == 2, "A must be a 2D array."

jacobien_e_pg = self.Get_jacobian_e_pg(matrixType)
weight_pg = self.Get_gauss(matrixType).weights
dN_e_pg = self.Get_dN_e_pg(matrixType)

DiffusePart_e_pg = np.einsum('ep,p,epki,kl,eplj->epij', jacobien_e_pg, weight_pg, dN_e_pg, A, dN_e_pg, optimize='optimal')
DiffusePart_e_pg = np.einsum('ep,p,epij->epji', jacobien_e_pg, weight_pg, dN_e_pg, optimize='optimal')

self.__dict_DiffusePart_e_pg[matrixType] = DiffusePart_e_pg

return self.__dict_DiffusePart_e_pg[matrixType].copy()

def Get_SourcePart_e_pg(self, matrixType: MatrixType) -> np.ndarray:
"""Returns the part that builds the source term (scalar).
SourcePart_e_pg = jacobian_e_pg, weight_pg, f_e_pg, N_pg'\n
SourcePart_e_pg = f_e_pg * jacobian_e_pg, weight_pg, N_pg'\n
Returns -> jacobian_e_pg, weight_pg, N_pg'
"""
Expand Down Expand Up @@ -1232,7 +1218,7 @@ def Get_Nodes_Line(self, line: Line) -> np.ndarray:
def Get_Nodes_Domain(self, domain: Domain) -> np.ndarray:
"""Returns nodes in the domain."""

assert isinstance(circle, Domain)
assert isinstance(domain, Domain)

coordo = self.coord

Expand Down Expand Up @@ -1442,7 +1428,7 @@ def Get_pointsInElem(self, coordinates: np.ndarray, elem: int) -> np.ndarray:
indexReord = np.append(np.arange(1, nPe), 0)
# Vectors i for edge segments
vect_i_b = coordo[connectMesh[indexReord]] - coordo[connectMesh]
# vect_i_b = np.einsum("ni,n->ni", vect_i_b, 1/np.linalg.norm(vect_i_b, axis=1), optimize="optimal")
vect_i_b = np.einsum("ni,n->ni", vect_i_b, 1/np.linalg.norm(vect_i_b, axis=1), optimize="optimal")

# normal vector to element face
vect_n = np.cross(vect_i_b[0], -vect_i_b[-1])
Expand Down Expand Up @@ -1564,20 +1550,19 @@ def __Get_Mapping(self, coordinates_n: np.ndarray, elements_e: np.ndarray, needC
# get groupElem datas
inDim = self.inDim
sysCoord_e = self.sysCoord_e # base change matrix for each element
matrixType = MatrixType.rigi
matrixType = MatrixType.mass
jacobian_e_pg = self.Get_jacobian_e_pg(matrixType, absoluteValues=False)
invF_e_pg = self.Get_invF_e_pg(matrixType)
dN_tild = self._dNtild()
nPg = invF_e_pg.shape[1]
gaussCoord_e_pg = self.Get_GaussCoordinates_e_p(matrixType)
xiOrigin = self.origin # origin of the reference element (xi, eta)

useIterative = False
# # Check whether iterative resolution is required
# # calculates the ratio between jacob min and max to detect if the element is distorted
# diff_e = jacobian_e_pg.max(1) * 1/jacobian_e_pg.min(1)
# error_e = 1 - diff_e # a perfect element has an error max <= 1e-12
# # a distorted element has a max error greater than 0
# useIterative = False
# Check whether iterative resolution is required
# calculates the ratio between jacob min and max to detect if the element is distorted
diff_e = jacobian_e_pg.max(1) * 1/jacobian_e_pg.min(1)
error_e = np.abs(1 - diff_e) # a perfect element has an error max <= 1e-12
# a distorted element has a max error greater than 0
useIterative_e = error_e > 1e-12
# useIterative = np.max(error_e) > 1e-12
else:
coordInElem_n = None
Expand Down Expand Up @@ -1622,36 +1607,27 @@ def ResearchFunction(e: int):

x0 = coordoElemBase[0,:dim] # orign of the real element (x,y)
xPs = coordinatesBase_n[:,:dim] # points coordinates (x,y)

# the fastest way, but may lead to shape functions that give values outside [0,1].
xiP: np.ndarray = xiOrigin + (xPs - x0) @ np.mean(invF_e_pg[e], 0)

# if not useIterative:
# if nPg == 1:
# # invF_e_pg is constant in the element
# xiP: np.ndarray = xiOrigin + (xPs - x0) @ invF_e_pg[e,0]
# else:
# # If the element have more than 1 integration point, it is necessary to choose the closest integration points. Because invF_e_pg is not constant in the element
# # for each node detected, we'll calculate its distance from all integration points and see where it's closest
# dist = np.zeros((xPs.shape[0], nPg))
# for p in range(nPg):
# dist[:,p] = np.linalg.norm(xPs - gaussCoord_e_pg[e, p, :dim], axis=1)
# invMin = invF_e_pg[e, np.argmin(dist, axis=1)]
# xiP: np.ndarray = xiOrigin + (xPs - x0) @ invMin
# else:
# # Here we need to construct the Jacobian matrices. This is the longest method here
# def Eval(xi: np.ndarray, xP):
# dN = _GroupElem._Evaluates_Functions(dN_tild, xi.reshape(1, -1))
# F = dN[0] @ coordoElemBase[:,:dim] # jacobian matrix
# J = x0 - xP + (xi - xiOrigin) @ F # cost function
# return J

# xiP = []
# for xP in xPs:
# res = least_squares(Eval, 0*xP, args=(xP,))
# xiP.append(res.x)

# xiP = np.array(xiP)

if not useIterative_e[e]:
# the fastest way
# available only for undistorted meshes !
# always valid for TRI3 and structured meshes.
xiP: np.ndarray = xiOrigin + (xPs - x0) @ invF_e_pg[e,0]

else:
# Here we need to construct the Jacobian matrices. This is the longest method here
def Eval(xi: np.ndarray, xP):
dN = _GroupElem._Evaluates_Functions(dN_tild, xi.reshape(1, -1))
F = dN[0] @ coordoElemBase[:,:dim] # jacobian matrix
J = x0 + (xi - xiOrigin) @ F - xP # cost function
return J

xiP = []
for xP in xPs:
res = least_squares(Eval, 0*xP, args=(xP,))
xiP.append(res.x)

xiP = np.array(xiP)

coordInElem_n[nodesInElement,:] = xiP.copy()

Expand Down
18 changes: 11 additions & 7 deletions EasyFEA/fem/_mesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,7 @@ def Translate(self, dx: float=0.0, dy: float=0.0, dz: float=0.0) -> None:
oldCoord = self.coordGlob
newCoord = oldCoord + np.array([dx, dy, dz])
for grp in self.dict_groupElem.values():
grp.coordGlob = newCoord
grp.coordGlob = newCoord
self._Notify('The mesh has been modified')


Expand Down Expand Up @@ -429,23 +429,23 @@ def Get_leftDispPart(self, matrixType: MatrixType) -> np.ndarray:

def Get_ReactionPart_e_pg(self, matrixType: MatrixType) -> np.ndarray:
"""Returns the part that builds the reaction term (scalar).
ReactionPart_e_pg = jacobian_e_pg * weight_pg * r_e_pg * N_pg' * N_pg\n
ReactionPart_e_pg = r_e_pg * jacobian_e_pg * weight_pg * N_pg' * N_pg\n
Returns -> jacobian_e_pg * weight_pg * N_pg' * N_pg
"""
return self.groupElem.Get_ReactionPart_e_pg(matrixType)

def Get_DiffusePart_e_pg(self, matrixType: MatrixType, A: np.ndarray) -> np.ndarray:
def Get_DiffusePart_e_pg(self, matrixType: MatrixType) -> np.ndarray:
"""Returns the part that builds the diffusion term (scalar).
DiffusePart_e_pg = jacobian_e_pg * weight_pg * k * dN_e_pg' * A * dN_e_pg\n
DiffusePart_e_pg = k_e_pg * jacobian_e_pg * weight_pg * dN_e_pg' * A * dN_e_pg\n
Returns -> jacobian_e_pg * weight_pg * dN_e_pg' * A * dN_e_pg
Returns -> jacobian_e_pg * weight_pg * dN_e_pg'
"""
return self.groupElem.Get_DiffusePart_e_pg(matrixType, A)
return self.groupElem.Get_DiffusePart_e_pg(matrixType)

def Get_SourcePart_e_pg(self, matrixType: MatrixType) -> np.ndarray:
"""Returns the part that builds the source term (scalar).
SourcePart_e_pg = jacobian_e_pg, weight_pg, f_e_pg, N_pg'\n
SourcePart_e_pg = f_e_pg * jacobian_e_pg, weight_pg, N_pg'\n
Returns -> jacobian_e_pg, weight_pg, N_pg'
"""
Expand Down Expand Up @@ -875,6 +875,10 @@ def Calc_projector(oldMesh: Mesh, newMesh: Mesh) -> sp.csr_matrix:
assert oldMesh.dim == newMesh.dim, "Mesh dimensions must be the same."

tic = Tic()

distoredMesh = np.max(np.abs(1 - oldMesh.Get_Quality('jacobian'))) > 1e-12
if distoredMesh:
Display.MyPrintError("Warning: distorted elements have been detected in the mesh.\nThey may lead to projection errors!")

detectedNodes, detectedElements_e, connect_e_n, coordo_n = oldMesh.groupElem.Get_Mapping(newMesh.coord)
# detectedNodes (size(connect_e_n)) are the nodes detected in detectedElements_e
Expand Down
5 changes: 3 additions & 2 deletions EasyFEA/materials/_elastic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -560,8 +560,9 @@ def _Behavior(self, dim: int=None, P: np.ndarray=None):

material_cM = Heterogeneous_Array(material_cM)

# # Verify that C = S^-1#
# assert np.linalg.norm(material_sM - np.linalg.inv(material_cM)) < 1e-10
# # Verify that S = C^-1
# assert np.linalg.norm(material_sM - np.linalg.inv(material_cM)) < 1e-10
# # Verify that C = S^-1
# assert np.linalg.norm(material_cM - np.linalg.inv(material_sM)) < 1e-10

# Performs a base change to orient the material in space
Expand Down
Loading

0 comments on commit 886852f

Please sign in to comment.