diff --git a/NonLinearSolver/internalPoints/ipIMDEACPUMAT.cpp b/NonLinearSolver/internalPoints/ipIMDEACPUMAT.cpp
index 82f24f03c1dfe00ef3393fc30156e710c83a148e..87be9731b9b6effdb13561a7174dfa1389c8decc 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 a051f44dfd50d14d9d13f8c54cfacf8686549854..2e5ad75af8370d8b0475ad77f54022d0d7f965fb 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 5a4d7a635994875a045ec071c7acf4409ef8059e..f86605eaa297fc2950659b120002b2ca8063ef0b 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 d33bd712b31eacd069ee85c0e8725bc49ec9c2bc..9924c44ef1fd4e5ed76ada7f19e26b3afb7b0a2c 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 2346dd18a392ce0da9945f908e4baf83756c2c0d..29dd09e4785554473cf4b46e81a25dd89a81626b 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 3193d7dbeb9c6312b81b19832925461c6b4631fa..0431aa83912ad90d910b28c192297b0f5a1dfc6e 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 5e7a920211841501622ee5093537fd04bd5f77f3..c032ef64e449d891e2cef6f97fb457c649efb801 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 0093669f87e1031bed91f1ceb9a9f26e3a5f9f31..785590dfaba180445b7b5157e1365e4d5757943d 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();