From 36dfe48715a5b8bd402c3929ec916a5f170fbe40 Mon Sep 17 00:00:00 2001
From: gabriel <gabriel.zarzoso@imdea.org>
Date: Fri, 23 Feb 2024 17:09:08 +0100
Subject: [PATCH] running with euler angles
---
.../internalPoints/ipIMDEACPUMAT.cpp | 6 ++--
NonLinearSolver/materialLaw/imdeacp.f | 4 ++-
.../materialLaw/library_crysplas_CAPSUL_v0.f | 24 +++++++-------
.../materialLaw/mlawIMDEACPUMAT.cpp | 12 +++++--
NonLinearSolver/materialLaw/mlawIMDEACPUMAT.h | 2 +-
dG3D/benchmarks/imdeaCP/cube.py | 32 +++++++++++++++----
dG3D/src/dG3DMaterialLaw.cpp | 4 +--
dG3D/src/dG3DMaterialLaw.h | 2 +-
8 files changed, 57 insertions(+), 29 deletions(-)
diff --git a/NonLinearSolver/internalPoints/ipIMDEACPUMAT.cpp b/NonLinearSolver/internalPoints/ipIMDEACPUMAT.cpp
index 82f24f03c..87be9731b 100644
--- a/NonLinearSolver/internalPoints/ipIMDEACPUMAT.cpp
+++ b/NonLinearSolver/internalPoints/ipIMDEACPUMAT.cpp
@@ -14,9 +14,9 @@
double IPIMDEACPUMAT::get(int comp) const
{
// for vizualization, with tag in ipField.h (only variables)
- //if(comp==IPField::PLASTICSTRAIN)
- // return _statev[18];
- //else
+ if(comp==IPField::PLASTICSTRAIN)
+ return _statev[45];
+ else
return IPUMATInterface::get(comp);
}
diff --git a/NonLinearSolver/materialLaw/imdeacp.f b/NonLinearSolver/materialLaw/imdeacp.f
index a051f44df..2e5ad75af 100755
--- a/NonLinearSolver/materialLaw/imdeacp.f
+++ b/NonLinearSolver/materialLaw/imdeacp.f
@@ -286,7 +286,7 @@ c READ STATE VARIABLES:
c write(*,*)'Start read state variables'
c If time=0 initialize variables
c write(*,*)'TIME',time(1),time(2)
-c write(*,*)'props',props(1)
+c write(*,*)'props',props(1),props(2),props(3),props(4)
if(time(2).EQ.0d0) THEN
@@ -966,6 +966,8 @@ C ROTATED_ORIENT=RR*ORIENT_MAT
C SAVE acumulated shear strain
index=index+4
statev(index)=gamma_tot_pred
+c write(*,*)'====',index,statev(index),'===='
+c write(*,*)'Def grad inc',DFGRD_INC
c write(*,*)'Save accumulated shear strain'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
diff --git a/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f b/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f
index 5a4d7a635..f86605eaa 100644
--- a/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f
+++ b/NonLinearSolver/materialLaw/library_crysplas_CAPSUL_v0.f
@@ -1554,12 +1554,12 @@ c write(*,*)' numprocess,',numprocess
real*8 orient_MAT(3,3),orient_mat_tr(3,3)
c write(*,*)'I am in lybrary_crysplas_CAPSUL cm3Libraries--------'
-c orient1(1)=props(1)
-c orient1(2)=props(2)
-c orient1(3)=props(3)
- orient1(1)=0.494651
- orient1(2)=0.869092
- orient1(3)=0.000000
+ orient1(1)=props(4)
+ orient1(2)=props(5)
+ orient1(3)=props(6)
+c orient1(1)=0.494651
+c orient1(2)=0.869092
+c orient1(3)=0.000000
c write(*,*)'orient',orient1
call norm(dnorm,orient1,3,1)
c write(*,*)'dnorm',dnorm
@@ -1567,12 +1567,12 @@ c write(*,*)'dnorm',dnorm
orient1(i)=orient1(i)/dnorm
enddo
c write(*,*)'orient after dnorm',orient1
-c orient2(1)=props(4)
-c orient2(2)=props(5)
-c orient2(3)=props(6)
- orient2(1)=-0.869092
- orient2(2)=0.494651
- orient2(3)=0.000000
+ orient2(1)=props(7)
+ orient2(2)=props(8)
+ orient2(3)=props(9)
+c orient2(1)=-0.869092
+c orient2(2)=0.494651
+c orient2(3)=0.000000
call norm(dnorm,orient2,3,1)
do,i=1,3
orient2(i)=orient2(i)/dnorm
diff --git a/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp
index d33bd712b..9924c44ef 100644
--- a/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp
+++ b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.cpp
@@ -29,7 +29,7 @@ extern "C" {
-mlawIMDEACPUMAT::mlawIMDEACPUMAT(const int num, const char *propName) : mlawUMATInterface(num,0.,propName)
+mlawIMDEACPUMAT::mlawIMDEACPUMAT(const int num, const char *propName, double euler0, double euler1, double euler2) : mlawUMATInterface(num,0.,propName)
{
//should depend on the umat: here we read the grains.inp file with 3 euler angle, 2 parameters equal to 1, and the number of internal variables (23)
std::ifstream in(propName);
@@ -49,8 +49,14 @@ mlawIMDEACPUMAT::mlawIMDEACPUMAT(const int num, const char *propName) : mlawUMA
i++;
}
-
-
+ props[i] = cos(euler0)*cos(euler2)-sin(euler0)*sin(euler2)*cos(euler1);
+ props[i+1] = sin(euler0)*cos(euler2)+cos(euler0)*sin(euler2)*cos(euler1);
+ props[i+2] = sin(euler2)*sin(euler1);
+ props[i+3] = -cos(euler0)*sin(euler2)-sin(euler0)*cos(euler2)*cos(euler1);
+ props[i+4] = -sin(euler0)*sin(euler2)+cos(euler0)*cos(euler2)*cos(euler1);
+ props[i+5] = cos(euler2)*sin(euler1);
+
+ std::cout << "Rotation Matrix -> umat :" << ","<<props[3] << ","<< props[4] << ","<< props[5] << ","<< props[6] << ","<< props[7] << ","<< props[8] << std::endl;
in.close();
// adapt propertie ratio and shear modulus and density
diff --git a/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.h b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.h
index 2346dd18a..29dd09e47 100644
--- a/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.h
+++ b/NonLinearSolver/materialLaw/mlawIMDEACPUMAT.h
@@ -72,7 +72,7 @@ class mlawIMDEACPUMAT : public mlawUMATInterface
{
public:
- mlawIMDEACPUMAT(const int num, const char *propName);
+ mlawIMDEACPUMAT(const int num, const char *propName, double euler0, double euler1, double euler2);
#ifndef SWIG
mlawIMDEACPUMAT(const mlawIMDEACPUMAT &source);
diff --git a/dG3D/benchmarks/imdeaCP/cube.py b/dG3D/benchmarks/imdeaCP/cube.py
index 3193d7dbe..0431aa839 100644
--- a/dG3D/benchmarks/imdeaCP/cube.py
+++ b/dG3D/benchmarks/imdeaCP/cube.py
@@ -1,17 +1,32 @@
#coding-Utf-8-*-
from gmshpy import *
-from dG3Dpy import*
+from dG3DpyDebug import*
from math import*
+import numpy as np
#script to launch PBC problem with a python script
-# material law
-lawnum = 11 # unique number of law
+angles_loaded=np.loadtxt('euler_angles.txt')
+print('angles_loaded 0',angles_loaded[0,:],'shape',angles_loaded.shape)
# material law
-lawnum = 11 # unique number of law
-law1 = IMDEACPUMATDG3DMaterialLaw(lawnum,'crystal.prop');
+
+#lawnum = 1 # unique number of law
+#law1 = IMDEACPUMATDG3DMaterialLaw(lawnum,'crystal.prop',angles_loaded[lawnum-1,0],angles_loaded[lawnum-1,1],angles_loaded[lawnum-1,2]);
+
+law=[]
+for i in range(len(angles_loaded[:,0])):
+ lawnum = i+1 # unique number of law
+ print('============= CRYSTAL ',lawnum ,'====================')
+ angle1=angles_loaded[lawnum-1,0]
+ angle2=angles_loaded[lawnum-1,1]
+ angle3=angles_loaded[lawnum-1,2]
+ law1 = IMDEACPUMATDG3DMaterialLaw(lawnum,'crystal.prop',angle1,angle2,angle3);
+ law.append(IMDEACPUMATDG3DMaterialLaw(lawnum,'crystal.prop',angle1,angle2,angle3))
+ print('law',str(i+1),'associated', 'angles:', angle1, angle2, angle3 )
+ print('----->law',i,law[i])
+ print('----->law1',i,law1)
#lawPert1 = dG3DMaterialLawWithTangentByPerturbation(law1, 1e-8)
@@ -42,7 +57,12 @@ nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
mysolver = nonLinearMechSolver(1000)
mysolver.loadModel(meshfile)
mysolver.addDomain(myfield1)
-mysolver.addMaterialLaw(law1)
+
+
+law2=law1
+mysolver.addMaterialLaw(law2)
+
+
#mysolver.addMaterialLaw(lawPert1)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 5e7a92021..c032ef64e 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -7766,10 +7766,10 @@ void VEVPUMATDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
//
//
-IMDEACPUMATDG3DMaterialLaw::IMDEACPUMATDG3DMaterialLaw(const int num, const char *propName) :
+IMDEACPUMATDG3DMaterialLaw::IMDEACPUMATDG3DMaterialLaw(const int num, const char *propName, double euler0, double euler1, double euler2) :
dG3DMaterialLaw(num,0.,true)
{
- _vevplaw = new mlawIMDEACPUMAT(num,propName);
+ _vevplaw = new mlawIMDEACPUMAT(num,propName,euler0,euler1,euler2);
double nu = _vevplaw->poissonRatio();
double mu = _vevplaw->shearModulus();
double E = 2.*mu*(1.+nu);
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 0093669f8..785590dfa 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -2609,7 +2609,7 @@ class IMDEACPUMATDG3DMaterialLaw : public dG3DMaterialLaw{
#endif //SWIG
public:
- IMDEACPUMATDG3DMaterialLaw(const int num, const char *propName);
+ IMDEACPUMATDG3DMaterialLaw(const int num, const char *propName, double euler0, double euler1, double euler2);
#ifndef SWIG
IMDEACPUMATDG3DMaterialLaw(const IMDEACPUMATDG3DMaterialLaw &source);
virtual ~IMDEACPUMATDG3DMaterialLaw();
--
GitLab