Correct fe2 constitutivelaw and update examples

examples/FE2/FE2.py

@@ -37,22 +37,20 @@
 fd.weakform.StressEquilibrium("FEM", name = "WeakForm") 
 
 #Create a global assembly
-fd.Assembly.create("WeakForm", "macro", name="Assembly", MeshChange = True) 
+macro_assembly = fd.Assembly.create("WeakForm", "macro", name="Assembly", MeshChange = True) 
 
 #Define a new static problem
 pb = fd.problem.NonLinear("Assembly")
 # Problem.set_nr_criterion("Displacement")
 
 #create a 'result' folder and set the desired ouputs
 if not(os.path.isdir('results')): os.mkdir('results')
-pb.add_output('results/FE2', 'Assembly', ['Disp', 'Stress', 'Strain', 'Stress_vm', 'Wm'], output_type='Node', file_format ='vtk')    
-pb.add_output('results/FE2', 'Assembly', ['Stress', 'Stress', 'Stress_vm'], output_type='Element', file_format ='vtk')    
+res = pb.add_output('results/FE2', 'Assembly', ['Disp', 'Stress', 'Strain', 'Wm'])    
 
 #output result for a random micro cell (here for the 5th integration point)
 #Warning : the results of micro_cells are automatically saved at each NR iteration (several iteration per time iteration)
-# Problem.initialize(0) #to build prolems
-# micro_cells.list_problem[5].add_output('results/micro_cell', micro_cells.list_assembly[5], ['disp', 'stress', 'strain', 'stress_vm'], output_type='Node', file_format ='vtk')    
-
+pb.initialize()
+res_micro = micro_cells.list_problem[5].add_output('results/micro_cell', micro_cells.list_assembly[5], ['Disp', 'Stress', 'Strain'])    
 
 #Definition of the set of nodes for boundary conditions
 crd = mesh_macro.nodes
@@ -68,12 +66,6 @@
 #displacement on right (ux=0.1mm)
 pb.bc.add('Dirichlet',right , 'DispX', 0.1) 
 
-pb.apply_boundary_conditions()
-
 #Solve problem
 pb.nlsolve()
 
-#---------- Post-Treatment ----------
-#Get the stress tensor, strain tensor, and displacement (nodal values)
-res_nd = pb.get_results("Assembly", ['Stress_VM','Strain'], 'Node')
-U = pb.get_disp()

examples/PeriodicBoundaryConditions/2D_Periodic_BC_in_composite_ply.py

@@ -1,5 +1,6 @@
 import fedoo as fd
 import numpy as np
+from fedoo.util.abaqus_inp import *
 import time
 
 #------------------------------------------------------------------------------
@@ -18,7 +19,7 @@
 # INP = Util.ReadINP('Job-1.inp')
 
 #warning: this mesh is not periodic and should be replaced by a better one !
-INP = fd.util.ReadINP('cell_taffetas.inp') 
+INP = ReadINP('cell_taffetas.inp') 
 mesh = INP.toMesh()
 
 E_fiber = 250e3
@@ -58,12 +59,12 @@
 #------------------------------------------------------------------------------
 #Mechanical weak formulation
 #------------------------------------------------------------------------------
-fd.weakform.StressEquilibrium("ElasticLaw")
+wf = fd.weakform.StressEquilibrium("ElasticLaw")
 
 #------------------------------------------------------------------------------
 #Global Matrix assembly
 #------------------------------------------------------------------------------
-assemb = fd.Assembly.create("ElasticLaw", mesh) 
+assemb = fd.Assembly.create(wf, mesh) 
 
 #------------------------------------------------------------------------------
 #Static problem based on the just defined assembly

examples/PeriodicBoundaryConditions/3D_periodic_BC.py

@@ -0,0 +1,87 @@
+import fedoo as fd
+import numpy as np
+
+#--------------- Pre-Treatment --------------------------------------------------------
+
+fd.ModelingSpace("3D")
+# mesh = fd.mesh.import_file('../../util/meshes/gyroid_per.vtk')
+mesh = fd.Mesh.read('../../util/meshes/gyroid_per.vtk')
+# assert (mesh.is_periodic())
+
+#Definition of the set of nodes for boundary conditions
+crd = mesh.nodes 
+xmax = np.max(crd[:,0]) ; xmin = np.min(crd[:,0])
+ymax = np.max(crd[:,1]) ; ymin = np.min(crd[:,1])
+zmax = np.max(crd[:,2]) ; zmin = np.min(crd[:,2])
+center = [np.linalg.norm(crd,axis=1).argmin()]
+
+StrainNodes = mesh.add_virtual_nodes(2) #add virtual nodes for macro strain
+
+#Material definition
+material = fd.constitutivelaw.ElasticIsotrop(1e5, 0.3)
+wf = fd.weakform.StressEquilibrium(material)
+
+#Assembly
+assembly = fd.Assembly.create(wf, mesh) 
+
+#Type of problem 
+pb = fd.problem.Linear(assembly)
+
+
+#Boundary conditions
+E = [0,0,0,0,0,1] #[EXX, EYY, EZZ, EXY, EXZ, EYZ]
+#StrainNode[0] - 'DispX' is a virtual dof for EXX
+#StrainNode[0] - 'DispY' is a virtual dof for EYY
+#StrainNode[0] - 'DispZ' is a virtual dof for EZZ        
+#StrainNode[1] - 'DispX' is a virtual dof for EXY        
+#StrainNode[1] - 'DispY' is a virtual dof for EXZ
+#StrainNode[1] - 'DispZ' is a virtual dof for EYZ
+
+#define periodic boundary condition
+list_strain_nodes = [StrainNodes[0], StrainNodes[0], StrainNodes[0],
+                      StrainNodes[1], StrainNodes[1], StrainNodes[1]]
+list_strain_var = ['DispX', 'DispY', 'DispZ','DispX', 'DispY', 'DispZ']
+
+bc_periodic = fd.constraint.PeriodicBC(list_strain_nodes, list_strain_var)
+pb.bc.add(bc_periodic)
+
+#boundary conditions
+pb.bc.add('Dirichlet',center,'DispX', 0)
+pb.bc.add('Dirichlet',center,'DispY', 0)
+pb.bc.add('Dirichlet',center,'DispZ', 0)
+
+pb.bc.add('Dirichlet',[StrainNodes[0]],'DispX', E[0]) #EpsXX
+pb.bc.add('Dirichlet',[StrainNodes[0]],'DispY', E[1]) #EpsYY
+pb.bc.add('Dirichlet',[StrainNodes[0]],'DispZ', E[2]) #EpsZZ
+pb.bc.add('Dirichlet',[StrainNodes[1]],'DispX', E[3]) #EpsXY
+pb.bc.add('Dirichlet',[StrainNodes[1]],'DispY', E[4]) #EpsXZ
+pb.bc.add('Dirichlet',[StrainNodes[1]],'DispZ', E[5]) #EpsYZ
+
+#lv--------------- Soe --------------------------------------------------------
+pb.solve()
+
+tensor_strain = assembly.sv['Strain']
+tensor_stress = assembly.sv['Stress']
+
+###############################################################################
+# print the macroscopic strain tensor and stress tensor
+mean_strain = [pb.get_disp('DispX')[-2], pb.get_disp('DispY')[-2], pb.get_disp('DispZ')[-2],
+               pb.get_disp('DispX')[-1], pb.get_disp('DispY')[-1], pb.get_disp('DispZ')[-1]]
+print('Strain tensor: ', 
+      np.array([[mean_strain[0], mean_strain[3], mean_strain[4]],
+                [mean_strain[3], mean_strain[1], mean_strain[5]],
+                [mean_strain[4], mean_strain[5], mean_strain[2]]])
+      )
+
+#Compute the mean stress tensor
+surf = mesh.bounding_box.volume #total surface of the domain = volume in 2d
+mean_stress = [1/surf*mesh.integrate_field(tensor_stress[i]) for i in range(6)]
+
+print('Stress tensor: ', 
+      np.array([[mean_stress[0], mean_stress[3], mean_stress[4]],
+                [mean_stress[3], mean_stress[1], mean_stress[5]],
+                [mean_stress[4], mean_stress[5], mean_stress[2]]])
+      )
+
+pb.get_results(assembly, 'Stress').plot('Stress', 'XX')
+