Skip to content
Snippets Groups Projects
Commit 98cf7fc3 authored by Van Dung NGUYEN's avatar Van Dung NGUYEN
Browse files

Merge branch 'vdg-cm3'

parents 3ccb4f1d 78a32517
No related branches found
No related tags found
No related merge requests found
...@@ -102,4 +102,5 @@ add_subdirectory(pathFollowing_TrueSnapback) ...@@ -102,4 +102,5 @@ add_subdirectory(pathFollowing_TrueSnapback)
add_subdirectory(pathFollowing_cohesive) add_subdirectory(pathFollowing_cohesive)
add_subdirectory(pathFollowingMultiscale) add_subdirectory(pathFollowingMultiscale)
add_subdirectory(multiscaleCohesiveTest2D_mixedMode) add_subdirectory(multiscaleCohesiveTest2D_mixedMode)
add_subdirectory(multiscaleCohesiveTest2D_fullDG)
add_subdirectory(compRVE) add_subdirectory(compRVE)
# test file
set(PYFILE model.py)
set(FILES2DELETE
E_*.msh
*.txt
*.csv
disp*
stress*
)
add_cm3python_mpi_test(4 ${PYFILE} "${FILES2DELETE}")
mm = 1.;
L = 1.*mm;
R = 0.15*L;
lsca1 = 0.2*L;
lsca2 = 0.4*lsca1;
Point(1) = {0,0,0,lsca1};
Point(2) = {-R,0,0,lsca2};
Point(3) = {-0.5*L,0,0,lsca2};
Point(4) = {-0.5*L,-0.5*L,0,lsca1};
Point(5) = {0,-0.5*L,0,lsca2};
Point(6) = {0,-R,0,lsca2};
Line(1) = {2, 3};
Line(2) = {3, 4};
Line(3) = {4, 5};
Line(4) = {5, 6};
Circle(5) = {2, 1, 6};
Line Loop(6) = {2, 3, 4, -5, 1};
Plane Surface(7) = {6};
Symmetry {1, 0, 0, 0} {
Duplicata { Surface{7}; }
}
Symmetry {0, 1, 0, 0} {
Duplicata { Surface{7, 8}; }
}
//Recombine Surface{14, 20, 8, 7};
Physical Line(1) = {3, 10};
Physical Line(2) = {9, 21};
Physical Line(3) = {22, 16};
Physical Line(4) = {15, 2};
Physical Surface(11) = {14, 20, 8, 7};
/*
Translate {0.1, 0, 0} {
Line{17, 24, 18, 5, 12, 4};
}
Translate {0, 0.1, 0} {
Line{18, 24, 13, 12, 5, 1};
}
*/
This diff is collapsed.
mm=1.0;
L=1.*mm;
H = 1*mm;
sl1=0.5*L;
Point(1)={0,0,0,sl1};
Point(2)={L,0,0,sl1};
Point(3)={L,H,0,sl1};
Point(4)={0,H,0,sl1};
Line(1)={1,2};
Line(2)={2,3};
Line(3)={3,4};
Line(4)={4,1};
l[0]=newreg;
Line Loop(l[0])={1,2,3,4};
Plane Surface(11)={l[0]};
Physical Surface(11)={11};
Physical Point(1)={1};
Physical Point(2)={2};
Physical Point(3)={3};
Physical Point(4)={4};
Transfinite Line {4, 2} = 2 Using Progression 1;
Transfinite Line {3, 1} = 2 Using Progression 1;
Transfinite Surface {11};
Recombine Surface {11};
Extrude {L, 0, 0} {
Line{2}; Layers{1}; Recombine;
}
Physical Surface(12) = {15};
Physical Line(1) = {13, 1};
Physical Line(2) = {12};
Physical Line(3) = {14, 3};
Physical Line(4) = {4};
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
6
1 0 0 0
2 1 0 0
3 1 1 0
4 0 1 0
5 2 0 0
6 2 1 0
$EndNodes
$Elements
12
1 15 2 1 1 1
2 15 2 2 2 2
3 15 2 3 3 3
4 15 2 4 4 4
5 1 2 1 1 1 2
6 1 2 3 3 3 4
7 1 2 4 4 4 1
8 1 2 2 12 5 6
9 1 2 1 13 2 5
10 1 2 3 14 3 6
11 3 2 11 11 1 2 3 4
12 3 2 12 15 2 3 6 5
$EndElements
#coding-Utf-8-*-
from gmshpy import *
from dG3Dpy import*
#script to launch PBC problem with a python script
#DEFINE MICRO PROBLEM
# micro-material law
lawnum1 = 11 # unique number of law
rho = 7850e-9
young = 28.9e3
nu = 0.3
sy0 = 99.
h = young/20.
harden = LinearExponentialJ2IsotropicHardening(1, sy0, h, 0., 10.)
cl = IsotropicCLengthLaw(1, 1e-2)
damlaw = SimpleSaturateDamageLaw(1, 50.,1.,0.,1.)
law1 = NonLocalDamageJ2HyperDG3DMaterialLaw(lawnum1,rho,young,nu,harden,cl,damlaw)
micromeshfile="micro.msh" # name of mesh file
#micromeshfile="square.msh" # name of mesh file
# creation of part Domain
nfield = 11 # number of the field (physical number of entity)
myfield1 = nonLocalDamageDG3DDomain(1000,nfield,0,lawnum1,0,1e3,2,1)
microBC = nonLinearPeriodicBC(1000,2)
microBC.setOrder(1)
microBC.setBCPhysical(1,2,3,4)
# periodiodic BC
method = 0 # Periodic mesh = 0 Langrange interpolation = 1 Cubic spline interpolation =2
degree = 5
addvertex = False
microBC.setPeriodicBCOptions(method, degree,addvertex)
# DEFINE MACROPROBLEM
matnum1 = 1;
macromat1 = dG3DMultiscaleMaterialLaw(matnum1, 1000)
macromat1.loadModel(micromeshfile);
macromat1.addDomain(myfield1)
macromat1.addMaterialLaw(law1);
macromat1.setPeriodicity(1.,0,0,"x")
macromat1.setPeriodicity(0,1.,0,"y")
macromat1.setPeriodicity(0,0,1.,"z")
macromat1.addMicroBC(microBC)
macromat1.setNumStep(3)
macromat1.setTolerance(1e-6,1e-10)
macromat1.setSystemType(1)
macromat1.Solver(2)
#for i in range(0,1):
macromat1.setViewAllMicroProblems(True,0)
macromat1.setBlockDamageAfterFailureOnset(True)
macromat1.Scheme(1)
macromat1.setSameStateCriterion(1e-16)
macromat1.stressAveragingFlag(True) # set stress averaging ON- 0 , OFF-1
macromat1.setStressAveragingMethod(0)
macromat1.tangentAveragingFlag(True) # set tangent averaging ON -0, OFF -1
macromat1.setTangentAveragingMethod(2,1e-6);
macromat1.internalPointBuildView("Green-Lagrange_xx",IPField.STRAIN_XX);
macromat1.internalPointBuildView("Green-Lagrange_yy",IPField.STRAIN_YY);
macromat1.internalPointBuildView("Green-Lagrange_xy",IPField.STRAIN_XY);
macromat1.internalPointBuildView("Green-Lagrange_zz",IPField.STRAIN_ZZ);
macromat1.internalPointBuildView("Green-Lagrange_yz",IPField.STRAIN_YZ);
macromat1.internalPointBuildView("Green-Lagrange_xz",IPField.STRAIN_XZ);
macromat1.internalPointBuildView("sig_xx",IPField.SIG_XX);
macromat1.internalPointBuildView("sig_yy",IPField.SIG_YY);
macromat1.internalPointBuildView("sig_xy",IPField.SIG_XY);
macromat1.internalPointBuildView("sig_xz",IPField.SIG_XZ);
macromat1.internalPointBuildView("sig_zz",IPField.SIG_ZZ);
macromat1.internalPointBuildView("sig_yz",IPField.SIG_YZ);
macromat1.internalPointBuildView("sig_VM",IPField.SVM);
macromat1.internalPointBuildView("F_xx",IPField.F_XX);
macromat1.internalPointBuildView("F_yy",IPField.F_YY);
macromat1.internalPointBuildView("F_xy",IPField.F_XY);
macromat1.internalPointBuildView("F_yx",IPField.F_YX);
macromat1.internalPointBuildView("P_xx",IPField.P_XX);
macromat1.internalPointBuildView("P_yy",IPField.P_YY);
macromat1.internalPointBuildView("P_xy",IPField.P_XY);
macromat1.internalPointBuildView("P_yx",IPField.P_YX);
macromat1.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN);
macromat1.internalPointBuildView("Damage",IPField.DAMAGE);
macromat1.dirichletBC("Face",11,2,0.)
lcohNum = 13
lawCoh = dG3DMultiscaleCohesiveLaw(lcohNum,0.7)
lawCoh.setCharacteristicLength(1.)
lawCoh.setExtractCohesiveLawFromMicroDamage(True)
lawCoh.setLostSolutionUniquenssTolerance(0.)
macromeshfile="model.msh" # name of mesh file
# solver
sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
soltype =1 # StaticLinear=0 (default) StaticNonLinear=1
nstep = 100 # number of step (used only if soltype=1)
ftime =1. # Final time (used only if soltype=1)
tol=1.e-5 # relative tolerance for NR scheme (used only if soltype=1)
nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
# creation of macro part Domain
dim =2
beta1 = 50;
fullDG = False;
averageStrainBased = True
# non DG domain
nfield1 = 11
macrodomain1 = dG3DDomain(10,nfield1,0,matnum1,fullDG,dim)
macrodomain1.stabilityParameters(beta1)
macrodomain1.setDistributedOtherRanks(True)
macrodomain1.distributeOnRootRank(False)
nfield2 = 12
macrodomain2 = dG3DDomain(11,nfield2,0,matnum1,fullDG,dim)
macrodomain2.stabilityParameters(beta1)
macrodomain2.setDistributedOtherRanks(True)
macrodomain2.distributeOnRootRank(False)
# interface domain
beta1 = 100;
interdomain1 = interDomainBetween3D(12,macrodomain1,macrodomain2,lcohNum,matnum1)
interdomain1.stabilityParameters(beta1)
interdomain1.averageStrainBased(averageStrainBased)
interdomain1.setDistributedOtherRanks(True)
interdomain1.distributeOnRootRank(False)
# creation of Solver
mysolver = nonLinearMechSolver(1000)
mysolver.loadModel(macromeshfile)
mysolver.addDomain(macrodomain1)
mysolver.addDomain(macrodomain2)
mysolver.addDomain(interdomain1)
mysolver.addMaterialLaw(macromat1)
mysolver.addMaterialLaw(lawCoh)
mysolver.setMultiscaleFlag(True)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol)
mysolver.options("-ksp_type preonly -pc_type lu -pc_factor_shift_type NONZERO -pc_factor_mat_solver_package petsc")
# boundary condition
mysolver.displacementBC("Face",11,2,0.0)
mysolver.displacementBC("Face",12,2,0.0)
mysolver.displacementBC("Face",11,1,0.0)
mysolver.displacementBC("Face",12,1,0.0)
mysolver.displacementBC("Edge",1,1,0.0)
mysolver.displacementBC("Edge",2,0,0.02)
mysolver.displacementBC("Edge",4,0,0.0)
# archivage
mysolver.internalPointBuildView("Green-Lagrange_xx",IPField.STRAIN_XX, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_yy",IPField.STRAIN_YY, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_zz",IPField.STRAIN_ZZ, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_xy",IPField.STRAIN_XY, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_yz",IPField.STRAIN_YZ, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange_xz",IPField.STRAIN_XZ, 1, 1);
mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1);
mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1);
mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1);
mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1);
mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1);
mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1);
mysolver.internalPointBuildView("sig_VM",IPField.SVM, 1, 1);
mysolver.internalPointBuildView("Green-Lagrange equivalent strain",IPField.GL_EQUIVALENT_STRAIN, 1, 1);
mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);
mysolver.internalPointBuildView("Damage",IPField.DAMAGE,1,1);
mysolver.archivingForceOnPhysicalGroup("Edge",4,0)
mysolver.archivingInterfaceElementIP(11,12,IPField.DISP_JUMP_X,4,0,1)
mysolver.archivingInterfaceElementIP(11,12,IPField.DISP_JUMP_Y,4,0,1)
mysolver.archivingInterfaceElementIP(11,12,IPField.COHESIVE_JUMP_X,4,0,1)
mysolver.archivingInterfaceElementIP(11,12,IPField.COHESIVE_JUMP_Y,4,0,1)
mysolver.archivingInterfaceElementIP(11,12,IPField.INCOMPATIBLE_STRAIN_X,4,0,1)
mysolver.archivingInterfaceElementIP(11,12,IPField.INCOMPATIBLE_STRAIN_Y,4,0,1)
# solve
mysolver.solve()
check = TestCheck()
check.equal(-8.808509e+00,mysolver.getArchivedForceOnPhysicalGroup("Edge", 4, 0),1.e-4)
...@@ -316,6 +316,300 @@ void GeneralMultiscaleBulkFollwedCohesive3DLaw::createIPVariable(IPVariable* &ip ...@@ -316,6 +316,300 @@ void GeneralMultiscaleBulkFollwedCohesive3DLaw::createIPVariable(IPVariable* &ip
} }
dG3DMultiscaleCohesiveLaw::dG3DMultiscaleCohesiveLaw(const int num, const double surfaceRatio):
GeneralMultiscaleBulkFollwedCohesive3DLaw(num,true),_surfaceReductionRatio(surfaceRatio){};
dG3DMultiscaleCohesiveLaw::dG3DMultiscaleCohesiveLaw(const dG3DMultiscaleCohesiveLaw& src):
GeneralMultiscaleBulkFollwedCohesive3DLaw(src),_surfaceReductionRatio(src._surfaceReductionRatio){}
bool dG3DMultiscaleCohesiveLaw::brokenCheck(IPVariable* ipv) const{
dG3DMultiscaleIPVariable* mipv = dynamic_cast<dG3DMultiscaleIPVariable*>(ipv);
if (mipv != NULL){
return mipv->solverIsBroken();
}
else{
return false;
}
};
void dG3DMultiscaleCohesiveLaw::checkCohesiveInsertion(IPStateBase* ipsm, IPStateBase* ipsp, const bool forcedInsert) const{
// check failure onset from both negative and postive IPStates
MultiscaleFractureCohesive3DIPVariable* fMinusCur = dynamic_cast<MultiscaleFractureCohesive3DIPVariable*>(ipsm->getState(IPStateBase::current));
const MultiscaleFractureCohesive3DIPVariable* fMinusPrev = dynamic_cast<const MultiscaleFractureCohesive3DIPVariable*>(ipsm->getState(IPStateBase::previous));
MultiscaleFractureCohesive3DIPVariable* fPlusCur = dynamic_cast<MultiscaleFractureCohesive3DIPVariable*>(ipsp->getState(IPStateBase::current));
const MultiscaleFractureCohesive3DIPVariable* fPlusPrev = dynamic_cast<const MultiscaleFractureCohesive3DIPVariable*>(ipsp->getState(IPStateBase::previous));
if (fMinusCur == NULL or fPlusCur == NULL or fMinusPrev==NULL or fPlusPrev==NULL){
Msg::Fatal("In MultiscaleCohesive3DLaw::initCohesiveLaw, MultiscaleFractureCohesive3DIPVariable is not used ");
}
if (fMinusPrev->isbroken() or fPlusPrev->isbroken()){
Msg::Warning("Previous IPs are broken");
fMinusCur->broken();
fPlusCur->broken();
return;
}
bool solverIsBrokenMinus = static_cast<const dG3DMultiscaleIPVariable*>(fMinusCur->getIPvBulk())->solverIsBroken();
bool solverIsBrokenPlus = static_cast<const dG3DMultiscaleIPVariable*>(fPlusCur->getIPvBulk())->solverIsBroken();
if (solverIsBrokenMinus or solverIsBrokenPlus){
printf("rank %d: cohesive element is inserted \n",Msg::GetCommRank());
fMinusCur->broken();
fPlusCur->broken();
// broken is checked via bulk terial law, see function stress
// initialize cohesive law in negative part
Cohesive3DIPVariable* coMinus = fMinusCur->getIPvFrac();
SVector3 zeroVec;
STensor3 zeroTen;
SVector3 normDir = fPlusCur->getConstRefToCurrentOutwardNormal();
normDir.normalize();
coMinus->initializeFracture(fMinusCur->getConstRefToJump(),0.,0.,0.,0.,0.,true,zeroTen,0.,normDir,zeroVec,
fMinusCur->getIPvBulk());
// positive part take a same value as negative cohesive part
Cohesive3DIPVariable* coPlus = fPlusCur->getIPvFrac();
coPlus->initializeFracture(fPlusCur->getConstRefToJump(),0.,0.,0.,0.,0.,true,zeroTen,0.,normDir,zeroVec,
fPlusCur->getIPvBulk());
}
else{
fMinusCur->nobroken();
fPlusCur->nobroken();
}
};
void dG3DMultiscaleCohesiveLaw::transferInterfaceDataToBulk(IPVariable* ipv, const IPVariable* ipvprev, const bool stiff) const{
FractureCohesive3DIPVariable *fipv = static_cast < FractureCohesive3DIPVariable *> (ipv);
const FractureCohesive3DIPVariable *fipvprev = static_cast < const FractureCohesive3DIPVariable *> (ipvprev);
if (fipvprev->isbroken()){
dG3DMultiscaleIPVariable* ipvMultiscale = dynamic_cast<dG3DMultiscaleIPVariable*>(fipv->getIPvBulk());
const dG3DMultiscaleIPVariable* ipvMultiscaleprev = dynamic_cast<const dG3DMultiscaleIPVariable*>(fipvprev->getIPvBulk());
MultiscaleCohesive3DIPVariable* cohIpv = dynamic_cast<MultiscaleCohesive3DIPVariable*>(fipv->getIPvFrac());
const MultiscaleCohesive3DIPVariable* cohIpvprev = dynamic_cast<const MultiscaleCohesive3DIPVariable*>(fipvprev->getIPvFrac());
nonLinearMechSolver* solver = NULL;
if (ipvMultiscale!= NULL){
solver = ipvMultiscale->getSolver();
}
else {
Msg::Fatal("solver is not available MultiscaleDGCohesive3DLaw::transferInterfaceDataToBulk");
}
const homogenizedData* homoData = solver->getHomogenizationState(IPStateBase::current);
const homogenizedData* homoDataPrev = solver->getHomogenizationState(IPStateBase::previous);
SVector3 refN = fipv->getConstRefToReferenceOutwardNormal();
refN.normalize();
const STensor3& Fb0 = fipvprev->getBulkDeformationGradient();
const STensor3& Fb = fipv->getBulkDeformationGradient();
const SVector3& ujump = fipv->getConstRefToJump();
const SVector3& ujump0 = fipvprev->getConstRefToJump();
double beta = homoData->getActiveDamageVolume()/(_surfaceReductionRatio*solver->getRVEVolume());
double lM = homoData->getAverageLocalizationBandWidth();
STensor3& F = fipv->getRefToDeformationGradient();
F = fipvprev->getConstRefToDeformationGradient();
SVector3 dFbN(0.);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
dFbN(i) += (Fb(i,j) - Fb0(i,j))*refN(j);
}
}
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
F(i,j) += Fb(i,j) - Fb0(i,j)-beta*dFbN(i)*refN(j) + (beta/lM)*(ujump(i) - ujump0(i))*refN(j);
}
}
if (stiff){
STensor43& dFinterfacedF = fipv->getRefToDInterfaceDeformationGradientDDeformationGradient();
STensor33& dFinterfacedJump = fipv->getRefToDInterfaceDeformationGradientDJump();
dFinterfacedF *= 0.;
dFinterfacedJump *= 0.;
static STensor3 I(1.);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
dFinterfacedJump(i,j,k) += I(i,k)*refN(j)*(beta/lM);
for (int l=0; l<3; l++){
dFinterfacedF(i,j,k,l) += I(i,k)*I(j,l) - beta*I(i,k)*refN(j)*refN(l);
}
}
}
}
}
solver->setLocalizationNormal(refN);
solver->setRVELengthInNormalDirection(_L,_surfaceReductionRatio);
}
};
void dG3DMultiscaleCohesiveLaw::stress(IPVariable* ipv, const IPVariable* ipvprev, const bool stiff, const bool checkfrac){
FractureCohesive3DIPVariable *fipv = static_cast < FractureCohesive3DIPVariable *> (ipv);
const FractureCohesive3DIPVariable *fipvprev = static_cast < const FractureCohesive3DIPVariable *> (ipvprev);
if (fipvprev->isbroken()){
dG3DMultiscaleIPVariable* mipv = static_cast<dG3DMultiscaleIPVariable*>(fipv->getIPvBulk());
const dG3DMultiscaleIPVariable* mipvprev = static_cast<const dG3DMultiscaleIPVariable*>(fipvprev->getIPvBulk());
MultiscaleCohesive3DIPVariable* cohmipv = dynamic_cast<MultiscaleCohesive3DIPVariable*>(fipv->getIPvFrac());
const MultiscaleCohesive3DIPVariable* cohmipvprev = dynamic_cast<const MultiscaleCohesive3DIPVariable*>(fipvprev->getIPvFrac());
nonLinearMechSolver* solver = mipv->getSolver();
const homogenizedData* homoData = solver->getHomogenizationState(IPStateBase::current);
const homogenizedData* homoDataPrev = solver->getHomogenizationState(IPStateBase::previous);
SVector3 refN = fipv->getConstRefToReferenceOutwardNormal();
refN.normalize();
const STensor3& P = fipv->getConstRefToFirstPiolaKirchhoffStress();
SVector3& T = fipv->getRefToInterfaceForce();
T *= 0.;
for (int i=0; i<3; i++) {
for (int j=0; j<3; j++){
T(i) += P(i,j)*refN(j);
}
}
SVector3& cjump = fipv->getRefToCohesiveJump();
cjump = homoData->getHomogenizedCohesiveJump();
double normalPart = dot(cjump,refN);
if (!cohmipvprev->ifTension()){
for (int i=0; i<3; i++)
cjump(i) -= normalPart*refN(i);
}
if (normalPart < 0){
cohmipv->setTensionFlag(false);
}
else{
cohmipv->setTensionFlag(true);
}
const double& lprev = homoDataPrev->getAverageLocalizationBandWidth();
double& irr = fipv->getRefToIrreversibleEnergy();
irr = homoData->getIrreversibleEnergy()*lprev;
if (stiff){
const STensor43& dFinterfacedF = fipv->getConstRefToDInterfaceDeformationGradientDDeformationGradient();
const STensor33& dFinterfacedJump = fipv->getConstRefToDInterfaceDeformationGradientDJump();
const STensor43& dPdFinf = homoData->getHomogenizedTangentOperator_F_F();
STensor33& dTdF = fipv->getRefToDInterfaceForceDDeformationGradient();
STensor3& dTdjump = fipv->getRefToDInterfaceForceDjump();
dTdF *= 0.;
dTdjump *= 0.;
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int p=0; p<3; p++){
for (int q=0; q<3; q++){
for (int k=0; k<3; k++){
dTdjump(i,k) += dPdFinf(i,j,p,q)*refN(j)*dFinterfacedJump(p,q,k);
for (int l=0; l<3; l++){
dTdF(i,k,l) += dPdFinf(i,j,p,q)*refN(j)*dFinterfacedF(p,q,k,l);
}
}
}
}
}
}
const STensor33& dcjumpdFinf = homoData->getHomogenizedTangentOperator_CohesiveJump_F();
STensor33& dcjumpDF = fipv->getRefToDCohesiveJumpDDeformationGradient();
STensor3& dcjumpDJump = fipv->getRefToDCohesiveJumpDJump();
STensor33 dcjumpdFinfNew(dcjumpdFinf);
if (!cohmipvprev->ifTension()){
STensor3 INN(1.);
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
INN(i,j) -= refN(i)*refN(j);
}
}
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
for (int l=0; l<3; l++){
dcjumpdFinfNew(i,j,k) += dcjumpdFinf(l,j,k)*INN(l,i);
}
}
}
}
}
dcjumpDF *= 0.;
dcjumpDJump*= 0.;
for (int i=0; i<3; i++){
for (int p=0; p<3; p++){
for (int q=0; q<3; q++){
for (int k=0; k<3; k++){
dcjumpDJump(i,k) += dcjumpdFinfNew(i,p,q)*dFinterfacedJump(p,q,k);
for (int l=0; l<3; l++){
dcjumpDF(i,k,l) += dcjumpdFinfNew(i,p,q)*dFinterfacedF(p,q,k,l);
}
}
}
}
}
if (homoData->getIrreversibleEnergyExtractionFlag()){
STensor3& dirrEnergDF = fipv->getRefToDIrreversibleEnergyDDeformationGradient();
SVector3& dirrEnergDJump = fipv->getRefToDIrreversibleEnergyDJump();
dirrEnergDF *=0.;
dirrEnergDJump*=0.;
const STensor3& dirrEnergdFinf = homoData->getHomogenizedTangentOperator_IrreversibleEnergy_F();
for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
for (int k=0; k<3; k++){
dirrEnergDJump(k) += dirrEnergdFinf(i,j)*dFinterfacedJump(i,j,k);
for (int l=0; l<3; l++){
dirrEnergDF(k,l) += dirrEnergdFinf(i,j)*dFinterfacedF(i,j,k,l);
}
}
}
}
}
}
}
};
TwoFieldMultiscaleCohesive3DLaw::TwoFieldMultiscaleCohesive3DLaw(const int num, const double surfaceRatio, const bool init): TwoFieldMultiscaleCohesive3DLaw::TwoFieldMultiscaleCohesive3DLaw(const int num, const double surfaceRatio, const bool init):
GeneralMultiscaleBulkFollwedCohesive3DLaw(num,init),_surfaceReductionRatio(surfaceRatio){}; GeneralMultiscaleBulkFollwedCohesive3DLaw(num,init),_surfaceReductionRatio(surfaceRatio){};
...@@ -349,7 +643,7 @@ void TwoFieldMultiscaleCohesive3DLaw::transferInterfaceDataToBulk(IPVariable* ip ...@@ -349,7 +643,7 @@ void TwoFieldMultiscaleCohesive3DLaw::transferInterfaceDataToBulk(IPVariable* ip
refN.normalize(); refN.normalize();
#if 1 #if 0
const SVector3& iJump0 = cohIpvprev->getConstRefToIncompatibleStrainAtFailureOnset(); const SVector3& iJump0 = cohIpvprev->getConstRefToIncompatibleStrainAtFailureOnset();
const SVector3& iJump = fipv->getConstRefToIncompatibleStrainVector(); const SVector3& iJump = fipv->getConstRefToIncompatibleStrainVector();
...@@ -390,8 +684,12 @@ void TwoFieldMultiscaleCohesive3DLaw::transferInterfaceDataToBulk(IPVariable* ip ...@@ -390,8 +684,12 @@ void TwoFieldMultiscaleCohesive3DLaw::transferInterfaceDataToBulk(IPVariable* ip
const STensor3& Fb0 = fipvprev->getBulkDeformationGradient(); const STensor3& Fb0 = fipvprev->getBulkDeformationGradient();
const STensor3& Fb = fipv->getBulkDeformationGradient(); const STensor3& Fb = fipv->getBulkDeformationGradient();
const SVector3& ujump = fipv->getConstRefToJump();
const SVector3& ujump0 = fipvprev->getConstRefToJump();
double beta = homoData->getActiveDamageVolume()/solver->getRVEVolume(); double beta = homoData->getActiveDamageVolume()/solver->getRVEVolume();
double lM = homoData->getAverageLocalizationBandWidth();
SVector3 dFbN(0.); SVector3 dFbN(0.);
for (int i=0; i<3; i++){ for (int i=0; i<3; i++){
...@@ -405,7 +703,7 @@ void TwoFieldMultiscaleCohesive3DLaw::transferInterfaceDataToBulk(IPVariable* ip ...@@ -405,7 +703,7 @@ void TwoFieldMultiscaleCohesive3DLaw::transferInterfaceDataToBulk(IPVariable* ip
for (int i=0; i<3; i++){ for (int i=0; i<3; i++){
for (int j=0; j<3; j++){ for (int j=0; j<3; j++){
F(i,j) += Fb(i,j) - Fb0(i,j)-beta*dFbN(i)*refN(j) + (iJump(i) - iJump0(i))*refN(j); F(i,j) += Fb(i,j) - Fb0(i,j)-beta*dFbN(i)*refN(j)+ beta*(ujump(i)-ujump0(i))*refN(j)/lM + (iJump(i) - iJump0(i))*refN(j);
//F(i,j) += (1.-beta)*(Fb(i,j) - Fb0(i,j)) + (iJump(i) - iJump0(i))*refN(j); //F(i,j) += (1.-beta)*(Fb(i,j) - Fb0(i,j)) + (iJump(i) - iJump0(i))*refN(j);
} }
} }
...@@ -423,6 +721,7 @@ void TwoFieldMultiscaleCohesive3DLaw::transferInterfaceDataToBulk(IPVariable* ip ...@@ -423,6 +721,7 @@ void TwoFieldMultiscaleCohesive3DLaw::transferInterfaceDataToBulk(IPVariable* ip
for (int j=0; j<3; j++){ for (int j=0; j<3; j++){
for (int k=0; k<3; k++){ for (int k=0; k<3; k++){
dFinterfacedInstrain(i,j,k) += I(i,k)*refN(j); dFinterfacedInstrain(i,j,k) += I(i,k)*refN(j);
dFinterfacedJump(i,j,k) += I(i,k)*refN(j)*(beta/lM);
for (int l=0; l<3; l++){ for (int l=0; l<3; l++){
//dFinterfacedF(i,j,k,l) += (1-beta)*I(i,k)*I(j,l); //dFinterfacedF(i,j,k,l) += (1-beta)*I(i,k)*I(j,l);
dFinterfacedF(i,j,k,l) += I(i,k)*I(j,l) - beta*I(i,k)*refN(j)*refN(l); dFinterfacedF(i,j,k,l) += I(i,k)*I(j,l) - beta*I(i,k)*refN(j)*refN(l);
......
...@@ -175,6 +175,22 @@ class GeneralMultiscaleBulkFollwedCohesive3DLaw : public GeneralBulkFollwedCohes ...@@ -175,6 +175,22 @@ class GeneralMultiscaleBulkFollwedCohesive3DLaw : public GeneralBulkFollwedCohes
#endif // SWIG #endif // SWIG
}; };
class dG3DMultiscaleCohesiveLaw : public GeneralMultiscaleBulkFollwedCohesive3DLaw {
protected:
double _surfaceReductionRatio;
public:
dG3DMultiscaleCohesiveLaw(const int num, const double surfaceRatio=1.);
#ifndef SWIG
dG3DMultiscaleCohesiveLaw(const dG3DMultiscaleCohesiveLaw& src);
virtual ~dG3DMultiscaleCohesiveLaw(){}
virtual bool brokenCheck(IPVariable* ipv) const;
virtual void stress(IPVariable* ipv, const IPVariable* ipvprev, const bool stiff=true, const bool checkfrac=true);
virtual void checkCohesiveInsertion(IPStateBase* ipsm, IPStateBase* ipsp, const bool forcedInsert= false) const;
virtual void transferInterfaceDataToBulk(IPVariable* ipv, const IPVariable* ipvprev, const bool stiff) const;
virtual materialLaw* clone() const {return new dG3DMultiscaleCohesiveLaw(*this);};
#endif // SWIG
};
class TwoFieldMultiscaleCohesive3DLaw : public GeneralMultiscaleBulkFollwedCohesive3DLaw{ class TwoFieldMultiscaleCohesive3DLaw : public GeneralMultiscaleBulkFollwedCohesive3DLaw{
protected: protected:
......
...@@ -656,10 +656,18 @@ void dG3DTwoFieldForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector< ...@@ -656,10 +656,18 @@ void dG3DTwoFieldForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<
normal.normalize(); normal.normalize();
const SVector3& ujump = ipvm->getConstRefToJump(); const SVector3& ujump = ipvm->getConstRefToJump();
// if fracture
bool broken = false;
const FractureCohesive3DIPVariable *ipvmfprev = dynamic_cast<const FractureCohesive3DIPVariable*>(ipvmprev);
if(ipvmfprev != NULL){
broken = ipvmfprev->isbroken();
}
if (!broken){
const STensor3& Pm = ipvm->getConstRefToFirstPiolaKirchhoffStress(); const STensor3& Pm = ipvm->getConstRefToFirstPiolaKirchhoffStress();
const STensor3& Pp = ipvp->getConstRefToFirstPiolaKirchhoffStress(); const STensor3& Pp = ipvp->getConstRefToFirstPiolaKirchhoffStress();
// consistent term // consistent term
SVector3 averagePN(0.); SVector3 averagePN(0.);
for (int j=0; j<3; j++){ for (int j=0; j<3; j++){
...@@ -680,14 +688,6 @@ void dG3DTwoFieldForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector< ...@@ -680,14 +688,6 @@ void dG3DTwoFieldForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<
} }
} }
// if fracture
bool broken = false;
const FractureCohesive3DIPVariable *ipvmfprev = dynamic_cast<const FractureCohesive3DIPVariable*>(ipvmprev);
if(ipvmfprev != NULL){
broken = ipvmfprev->isbroken();
}
if (!broken){
// DG formulation // DG formulation
const STensor43& eHm = ipvm->getConstRefToElasticTangentModuli(); const STensor43& eHm = ipvm->getConstRefToElasticTangentModuli();
const STensor43& eHp = ipvp->getConstRefToElasticTangentModuli(); const STensor43& eHp = ipvp->getConstRefToElasticTangentModuli();
...@@ -762,6 +762,23 @@ void dG3DTwoFieldForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector< ...@@ -762,6 +762,23 @@ void dG3DTwoFieldForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<
const FractureCohesive3DIPVariable *ipvmf = static_cast<const FractureCohesive3DIPVariable*>(ipvm); const FractureCohesive3DIPVariable *ipvmf = static_cast<const FractureCohesive3DIPVariable*>(ipvm);
const FractureCohesive3DIPVariable *ipvpf = static_cast<const FractureCohesive3DIPVariable*>(ipvp); const FractureCohesive3DIPVariable *ipvpf = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
SVector3 T = ipvmf->getConstRefToInterfaceForce();
T += ipvpf->getConstRefToInterfaceForce();
T*= 0.5;
// Assembly consistency term
for(int j=0;j<nbFFm;j++){
for(int k=0;k<3;k++){
m(j+k*nbFFm) -= T(k)*(Valsm[j+0*nbFFm]*wJ);
}
}
for(int j=0;j<nbFFp;j++){
for(int k=0;k<3;k++){
m(j+k*nbFFp+nbdofm) += T(k)*(Valsp[j+0*nbFFp]*wJ);
}
}
SVector3 cjump = ipvmf->getConstRefToCohesiveJump(); SVector3 cjump = ipvmf->getConstRefToCohesiveJump();
cjump += ipvpf->getConstRefToCohesiveJump(); cjump += ipvpf->getConstRefToCohesiveJump();
cjump *= (0.5); cjump *= (0.5);
...@@ -844,6 +861,16 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri ...@@ -844,6 +861,16 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri
normal.normalize(); normal.normalize();
const SVector3& ujump = ipvm->getConstRefToJump(); const SVector3& ujump = ipvm->getConstRefToJump();
// if fracture
bool broken = false;
const FractureCohesive3DIPVariable *ipvmfprev = dynamic_cast<const FractureCohesive3DIPVariable*>(ipvmprev);
if(ipvmfprev != NULL){
broken = ipvmfprev->isbroken();
}
if (!broken){
const STensor43& Hm = ipvm->getConstRefToTangentModuli(); const STensor43& Hm = ipvm->getConstRefToTangentModuli();
const STensor43& Hp = ipvp->getConstRefToTangentModuli(); const STensor43& Hp = ipvp->getConstRefToTangentModuli();
...@@ -896,14 +923,6 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri ...@@ -896,14 +923,6 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri
} }
} }
// if fracture
bool broken = false;
const FractureCohesive3DIPVariable *ipvmfprev = dynamic_cast<const FractureCohesive3DIPVariable*>(ipvmprev);
if(ipvmfprev != NULL){
broken = ipvmfprev->isbroken();
}
if (!broken){
// DG formulation // DG formulation
const STensor43& eHm = ipvm->getConstRefToElasticTangentModuli(); const STensor43& eHm = ipvm->getConstRefToElasticTangentModuli();
const STensor43& eHp = ipvp->getConstRefToElasticTangentModuli(); const STensor43& eHp = ipvp->getConstRefToElasticTangentModuli();
...@@ -1017,24 +1036,97 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri ...@@ -1017,24 +1036,97 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri
} }
else{ else{
// twofield formulation const FractureCohesive3DIPVariable *ipvmf = static_cast<const FractureCohesive3DIPVariable*>(ipvm);
STensor3 NDaveragePNDinStrain(0.); const FractureCohesive3DIPVariable *ipvpf = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
for (int j=0; j<3; j++){
const STensor33& dTdFm = ipvmf->getConstRefToDInterfaceForceDDeformationGradient();
const STensor33& dTdFp = ipvpf->getConstRefToDInterfaceForceDDeformationGradient();
for(int j=0;j<nbFFm;j++){
for(int k=0;k<3;k++){ for(int k=0;k<3;k++){
for (int p=0; p<3; p++){ for (int p=0; p<nbFFm; p++){
for (int q=0; q< 3; q++){
for (int r=0; r<3; r++){
m(j+k*nbFFm, p+q*nbFFm) -= 0.5*dTdFm(k,q,r)*Gradsm[p+q*nbFFm](r)*(Valsm[j+k*nbFFm]*wJ);
}
}
}
for (int p=0; p<nbFFp; p++){
for (int q=0; q< 3; q++){
for (int r=0; r<3; r++){
m(j+k*nbFFm, p+q*nbFFp+nbdofm) -= 0.5*dTdFp(k,q,r)*Gradsp[p+q*nbFFp](r)*(Valsm[j+k*nbFFm]*wJ);
}
}
}
}
}
for(int j=0;j<nbFFp;j++){
for(int k=0;k<3;k++){
for (int p=0; p<nbFFm; p++){
for (int q=0; q< 3; q++){
for (int r=0; r<3; r++){
m(j+k*nbFFp+nbdofm,p+q*nbFFm) += 0.5*dTdFm(k,q,r)*Gradsm[p+q*nbFFm](r)*(Valsp[j+k*nbFFp]*wJ);
}
}
}
for (int p=0; p<nbFFp; p++){
for (int q=0; q< 3; q++){
for (int r=0; r<3; r++){
m(j+k*nbFFp+nbdofm,p+q*nbFFp+nbdofm) += 0.5*dTdFp(k,q,r)*Gradsp[p+q*nbFFp](r)*(Valsp[j+k*nbFFp]*wJ);
}
}
}
}
}
STensor3 dTDjump(ipvmf->getConstRefToDInterfaceForceDjump());
dTDjump += ipvpf->getConstRefToDInterfaceForceDjump();
dTDjump *= 0.5;
for(int j=0;j<nbFFm;j++){
for(int k=0;k<3;k++){
for (int p=0; p< nbFFm; p++){
for (int q=0; q<3; q++){ for (int q=0; q<3; q++){
NDaveragePNDinStrain(j,p) += 0.5*(Hm(j,k,p,q)+Hp(j,k,p,q))*normal(k)*normal(q); m(j+k*nbFFm,p+q*nbFFm) += dTDjump(k,q)*Valsm[p+0*nbFFm]*(Valsm[j+0*nbFFm]*wJ);
} }
}
for (int p=0; p< nbFFp; p++){
for (int q=0; q<3; q++){
m(j+k*nbFFm,p+q*nbFFm+nbdofm) -= dTDjump(k,q)*Valsp[p+0*nbFFp]*(Valsm[j+0*nbFFm]*wJ);
} }
} }
} }
}
for(int j=0;j<nbFFp;j++){
for(int k=0;k<3;k++){
for (int p=0; p<nbFFm; p++){
for (int q=0; q< 3; q++){
m(j+k*nbFFp+nbdofm, p+q*nbFFm) -= dTDjump(k,q)*Valsm[p+0*nbFFm]*(Valsp[j+0*nbFFp]*wJ);
}
}
for (int p=0; p<nbFFp; p++){
for (int q=0; q< 3; q++){
m(j+k*nbFFp+nbdofm, p+q*nbFFp+nbdofm) += dTDjump(k,q)*Valsp[p+0*nbFFp]*(Valsp[j+0*nbFFp]*wJ);
}
}
}
}
// twofield formulation
STensor3 dTdinStrain(ipvmf->getConstRefToDInterfaceForceDIncompatibleStrainVector());
dTdinStrain += ipvpf->getConstRefToDInterfaceForceDIncompatibleStrainVector();
dTdinStrain *= (0.5);
if (_dom->getEnrichedUnknownLocation() == partDomain::NODE){ if (_dom->getEnrichedUnknownLocation() == partDomain::NODE){
for(int j=0;j<nbFFm;j++){ for(int j=0;j<nbFFm;j++){
for(int k=0;k<3;k++){ for(int k=0;k<3;k++){
for (int p=0; p< nbFF; p++){ for (int p=0; p< nbFF; p++){
for (int q=0; q< 3; q++){ for (int q=0; q< 3; q++){
m(j+k*nbFFm,p+q*nbFF+nbdofm+nbdofp) -= NDaveragePNDinStrain(k,q)*Vals[p+q*nbFF]*(Valsm[j+k*nbFFm]*wJ); m(j+k*nbFFm,p+q*nbFF+nbdofm+nbdofp) -= dTdinStrain(k,q)*Vals[p+q*nbFF]*(Valsm[j+k*nbFFm]*wJ);
} }
} }
} }
...@@ -1043,7 +1135,7 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri ...@@ -1043,7 +1135,7 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri
for(int k=0;k<3;k++){ for(int k=0;k<3;k++){
for (int p=0; p< nbFF; p++){ for (int p=0; p< nbFF; p++){
for (int q=0; q< 3; q++){ for (int q=0; q< 3; q++){
m(j+k*nbFFp+nbdofm, p+q*nbFF+nbdofm+nbdofp) += NDaveragePNDinStrain(k,q)*Vals[p+q*nbFF]*(Valsp[j+k*nbFFp]*wJ); m(j+k*nbFFp+nbdofm, p+q*nbFF+nbdofm+nbdofp) += dTdinStrain(k,q)*Vals[p+q*nbFF]*(Valsp[j+k*nbFFp]*wJ);
} }
} }
} }
...@@ -1053,7 +1145,7 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri ...@@ -1053,7 +1145,7 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri
for(int j=0;j<nbFFm;j++){ for(int j=0;j<nbFFm;j++){
for(int k=0;k<3;k++){ for(int k=0;k<3;k++){
for (int q=0; q< 3; q++){ for (int q=0; q< 3; q++){
m(j+k*nbFFm,i+q*npts+nbdofm+nbdofp) -= NDaveragePNDinStrain(k,q)*(Valsm[j+k*nbFFm]*wJ); m(j+k*nbFFm,i+q*npts+nbdofm+nbdofp) -= dTdinStrain(k,q)*(Valsm[j+k*nbFFm]*wJ);
} }
} }
...@@ -1061,21 +1153,22 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri ...@@ -1061,21 +1153,22 @@ void dG3DTwoFieldStiffnessInter::get(MElement *ele,int npts,IntPt *GP, fullMatri
for(int j=0;j<nbFFp;j++){ for(int j=0;j<nbFFp;j++){
for(int k=0;k<3;k++){ for(int k=0;k<3;k++){
for (int q=0; q< 3; q++){ for (int q=0; q< 3; q++){
m(j+k*nbFFp+nbdofm, i+q*npts+nbdofm+nbdofp) += NDaveragePNDinStrain(k,q)*(Valsp[j+k*nbFFp]*wJ); m(j+k*nbFFp+nbdofm, i+q*npts+nbdofm+nbdofp) += dTdinStrain(k,q)*(Valsp[j+k*nbFFp]*wJ);
} }
} }
} }
} }
const FractureCohesive3DIPVariable *ipvmf = static_cast<const FractureCohesive3DIPVariable*>(ipvm);
const FractureCohesive3DIPVariable *ipvpf = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
const STensor33& SFm = ipvmf->getConstRefToDCohesiveJumpDDeformationGradient(); const STensor33& SFm = ipvmf->getConstRefToDCohesiveJumpDDeformationGradient();
const STensor33& SFp = ipvpf->getConstRefToDCohesiveJumpDDeformationGradient(); const STensor33& SFp = ipvpf->getConstRefToDCohesiveJumpDDeformationGradient();
const STensor3& Sum = ipvmf->getConstRefToDCohesiveJumpDJump();
const STensor3& Sup = ipvpf->getConstRefToDCohesiveJumpDJump();
STensor3 DdiffJumpDjump(1.); STensor3 DdiffJumpDjump(1.);
DdiffJumpDjump.daxpy(Sum,-0.5);
DdiffJumpDjump.daxpy(Sup,-0.5);
STensor3 DdiffJumpDinStrain = ipvmf->getConstRefToDCohesiveJumpDIncompatibleStrainVector(); STensor3 DdiffJumpDinStrain = ipvmf->getConstRefToDCohesiveJumpDIncompatibleStrainVector();
DdiffJumpDinStrain += ipvpf->getConstRefToDCohesiveJumpDIncompatibleStrainVector(); DdiffJumpDinStrain += ipvpf->getConstRefToDCohesiveJumpDIncompatibleStrainVector();
DdiffJumpDinStrain *= (-0.5); DdiffJumpDinStrain *= (-0.5);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment