diff --git a/dG3D/src/dG3DTerms.cpp b/dG3D/src/dG3DTerms.cpp
index 870fa8ac1beacabd8b7ec463ff238d318f761bf0..c24253e0108a60fc4eb9a3b35584739b1036b89e 100644
--- a/dG3D/src/dG3DTerms.cpp
+++ b/dG3D/src/dG3DTerms.cpp
@@ -796,6 +796,7 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
// Initialization of some data
int nbdof = g3DBilinearTerm<double,double>::space1.getNumKeys(ele);
int nbFF = ele->getNumShapeFunctions();
+ int numEdgeShapeFunc = 1*ele->getNumEdges();
const dG3DMaterialLaw *mlawbulk;
if(_mlaw->getType() == materialLaw::fracture)
mlawbulk = static_cast<const dG3DMaterialLaw* >((static_cast<const FractureByCohesive3DLaw* > (_mlaw) )->getBulkLaw());
@@ -1047,6 +1048,36 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
}
}
}
+
+ // MECA-CURL DATA
+ if (getNumConstitutiveCurlVariable() > 0)
+ {
+ if( ipv->getNumConstitutiveCurlVariable() > getNumConstitutiveCurlVariable() )
+ Msg::Error("Your material law uses more constitutive curl data variables than your domain dG3DStiffnessBulk::get");
+
+ // Basis functions for vector field at edges
+ const std::vector<SVector3> &curlVals = ipv->fcurl(&space1,ele,GP[i],"HcurlLegendre");
+ // Basis functions for Curl of vector field at edges
+ const std::vector<SVector3> &CurlcurlVals = ipv->Curlfcurl(&space1,ele,GP[i],"CurlHcurlLegendre");
+
+ for (unsigned int curlField = 0; curlField < ipv->getNumConstitutiveCurlVariable(); ++curlField)
+ {
+ // dPdB
+ const STensor33 & dPdVectorCurl = ipv->getConstRefTodPdVectorCurl()[curlField];
+ for (unsigned int j = 0; j < nbFF; ++j)
+ {
+ for (unsigned int k = 0; k < numEdgeShapeFunc; ++k)
+ {
+ const unsigned int indexCurl = 3+getNumNonLocalVariable()+getNumConstitutiveExtraDofDiffusionVariable()+curlField;
+ for (unsigned int kk = 0; kk < 3; ++kk)
+ for (unsigned int m = 0; m < 3; ++m)
+ {
+ //mStiff(j+kk*nbFF,k+indexCurl*numEdgeShapeFunc) += ratio * 0.0;
+ }
+ }
+ }
+ }
+ }
}
//
@@ -1200,6 +1231,31 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
}
}
}
+
+ // NON-LOCAL - CURL FIELD
+ if (getNumConstitutiveCurlVariable() > 0)
+ {
+ if( ipv->getNumConstitutiveCurlVariable() > getNumConstitutiveCurlVariable() )
+ Msg::Error("Your material law uses more constitutive curl data variables than your domain dG3DStiffnessBulk::get");
+
+ // Basis functions for vector field at edges
+ const std::vector<SVector3> &curlVals = ipv->fcurl(&space1,ele,GP[i],"HcurlLegendre");
+ // Basis functions for Curl of vector field at edges
+ const std::vector<SVector3> &CurlcurlVals = ipv->Curlfcurl(&space1,ele,GP[i],"CurlHcurlLegendre");
+
+ for (unsigned int curlField = 0; curlField < ipv->getNumConstitutiveCurlVariable(); ++curlField)
+ {
+ const SVector3 & dLocalVariabledVectorCurl = ipv->getConstRefTodLocalVariabledVectorCurl()[curlField];
+ for (unsigned int j = 0; j < nbFF; ++j)
+ {
+ for (unsigned int k = 0; k < numEdgeShapeFunc; ++k)
+ {
+ const unsigned int IndexCurl = 3+getNumNonLocalVariable()+getNumConstitutiveExtraDofDiffusionVariable()+curlField;
+ //mStiff() += ratio * 0.0;
+ }
+ }
+ }
+ }
}
}
}
@@ -1452,6 +1508,45 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
}
}
}
+
+ // EXTRA-DOF - CURL FIELD
+ if (getNumConstitutiveCurlVariable() > 0)
+ {
+ if( ipv->getNumConstitutiveCurlVariable() > getNumConstitutiveCurlVariable() )
+ Msg::Error("Your material law uses more constitutive curl data variables than your domain dG3DStiffnessBulk::get");
+
+ // Basis functions for vector field at edges
+ const std::vector<SVector3> &curlVals = ipv->fcurl(&space1,ele,GP[i],"HcurlLegendre");
+ // Basis functions for Curl of vector field at edges
+ const std::vector<SVector3> &CurlcurlVals = ipv->Curlfcurl(&space1,ele,GP[i],"CurlHcurlLegendre");
+
+ for (unsigned int extraDOFField = 0; extraDOFField < ipv->getNumConstitutiveExtraDofDiffusionVariable(); ++extraDOFField)
+ for (unsigned int curlField = 0; curlField < ipv->getNumConstitutiveCurlVariable(); ++curlField)
+ {
+ const STensor3 &dFluxdVectorCurl = ipv->getConstRefTodFluxdVectorCurl()[extraDOFField][curlField];
+ SVector3 dFieldSourcedVectorCurl(0.0);
+ SVector3 dMechSourcedVectorCurl(0.0);
+
+ if (getConstitutiveExtraDofDiffusionAccountFieldSource())
+ {
+ dFieldSourcedVectorCurl += ipv->getConstRefTodFieldSourcedVectorCurl()[extraDOFField][curlField];
+ }
+
+ if (getConstitutiveExtraDofDiffusionAccountMecaSource())
+ {
+ dMechSourcedVectorCurl += ipv->getConstRefTodMechanicalSourcedVectorCurl()[extraDOFField][curlField];
+ }
+
+ for (unsigned int j = 0; j < nbFF; ++j)
+ for (unsigned int k = 0; k < numEdgeShapeFunc; ++k)
+ {
+ const unsigned int IndexField = 3+getNumNonLocalVariable()+extraDOFField;
+ const unsigned int IndexCurl = 3+getNumNonLocalVariable()+getNumConstitutiveExtraDofDiffusionVariable()+curlField;
+
+ //mStiff;
+ }
+ }
+ }
}
}
@@ -1461,14 +1556,122 @@ void dG3DStiffnessBulk::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
* */
if (getNumConstitutiveCurlVariable() >0)
{
- // make curl-curl term
-
- // make curl-disp term
-
- // make curl-nonlocal term
-
+ for (int i = 0; i < npts; i++)
+ {
+ const IPStateBase *ips = (*vips)[i];
+ const dG3DIPVariableBase *ipv = static_cast<const dG3DIPVariableBase*>(ips->getState(IPStateBase::current));
+
+ if( ipv->getNumConstitutiveCurlVariable() > getNumConstitutiveCurlVariable() )
+ Msg::Error("Your material law uses more constitutive curl data variables than your domain dG3DStiffnessBulk::get");
+
+ const double weight = GP[i].weight;
+ const double detJ = ipv->getJacobianDeterminant(ele,GP[i]);
+ const double ratio = detJ*weight;
- // make curl- extraDof term
+ // Basis functions for vector field at edges
+ const std::vector<SVector3> &curlVals = ipv->fcurl(&space1,ele,GP[i],"HcurlLegendre");
+ // Basis functions for Curl of vector field at edges
+ const std::vector<SVector3> &CurlcurlVals = ipv->Curlfcurl(&space1,ele,GP[i],"CurlHcurlLegendre");
+
+ for (unsigned int curlField = 0; curlField < ipv->getNumConstitutiveCurlVariable(); ++curlField)
+ {
+ // make curl-curl term
+ for (unsigned int curlField2 = 0; curlField2 < ipv->getNumConstitutiveCurlVariable(); ++curlField2)
+ {
+ const STensor3 & dHdB = ipv->getConstRefTodVectorFielddVectorCurl()[curlField][curlField2];
+
+ STensor3 dSourceVectorFielddB(0.0);
+ STensor3 dMechFieldSourcedB(0.0);
+
+ if (getConstitutiveCurlAccountFieldSource())
+ {
+ dSourceVectorFielddB += ipv->getConstRefTodSourceVectorFielddVectorCurl()[curlField][curlField2];
+ }
+ if (getConstitutiveCurlAccountMecaSource())
+ {
+ dMechFieldSourcedB += ipv->getConstRefTodMechanicalFieldSourcedVectorCurl()[curlField][curlField2];
+ }
+
+ for (unsigned int j = 0; j < numEdgeShapeFunc; ++j)
+ for (unsigned int k = 0; k < numEdgeShapeFunc; ++k)
+ {
+ const unsigned int IndexCurl = 3+getNumNonLocalVariable()+getNumConstitutiveExtraDofDiffusionVariable()+curlField;
+ const unsigned int IndexCurl2 = 3+getNumNonLocalVariable()+getNumConstitutiveExtraDofDiffusionVariable()+curlField2;
+
+ //mStiff;
+ }
+ }
+
+ // make curl-disp term
+
+ // make curl-nonlocal term
+ if (getNumNonLocalVariable() > 0)
+ {
+ if( ipv->getNumberNonLocalVariable() > getNumNonLocalVariable() )
+ Msg::Error("Your material law uses more non local variables than your domain dG3DStiffnessBulk::get");
+
+ for (unsigned int nlk = 0; nlk < ipv->getNumberNonLocalVariable(); ++nlk)
+ {
+ const SVector3 & dVectorFielddNonLocalVar = ipv->getConstRefTodVectorFielddNonLocalVariable()[curlField][nlk];
+
+ SVector3 dSourceVectorFielddNonLocalVar(0.0);
+ SVector3 dMechFieldSourcedNonLocalVar(0.0);
+
+ if (getConstitutiveCurlAccountFieldSource())
+ {
+ dSourceVectorFielddNonLocalVar += ipv->getConstRefTodSourceVectorFielddNonLocalVariable()[curlField][nlk];
+ }
+ if (getConstitutiveCurlAccountMecaSource())
+ {
+ dMechFieldSourcedNonLocalVar += ipv->getConstRefTodMechanicalFieldSourcedNonLocalVariable()[curlField][nlk];
+ }
+
+ for (unsigned int j = 0; j < numEdgeShapeFunc; ++j)
+ for (unsigned int k = 0; k < nbFF; ++k)
+ {
+ const unsigned int IndexCurl = 3+getNumNonLocalVariable()+getNumConstitutiveExtraDofDiffusionVariable()+curlField;
+ const unsigned int IndexNlk = 3+nlk;
+
+ //mStiff;
+ }
+ }
+ }
+
+ // make curl- extraDof term
+ if (getNumConstitutiveExtraDofDiffusionVariable() > 0)
+ {
+ if (ipv->getNumConstitutiveExtraDofDiffusionVariable() > getNumConstitutiveExtraDofDiffusionVariable())
+ Msg::Error("Your material law uses more constitutive extra dof variables than your domain dG3DStiffnessBulk::get");
+
+ for (unsigned int extraDOFField = 0; extraDOFField < ipv->getNumConstitutiveExtraDofDiffusionVariable(); ++extraDOFField)
+ {
+ const SVector3 & dVectorFielddExtraDofField = ipv->getConstRefTodVectorFielddExtraDofField()[curlField][extraDOFField];
+ const STensor3 & dVectorFielddGradExtraDofField = ipv->getConstRefTodVectorFielddGradExtraDofField()[curlField][extraDOFField];
+
+ SVector3 dSourceVectorFielddExtraDofField(0.0);
+ SVector3 dMechFieldSourcedExtraDofField(0.0);
+
+ if (getConstitutiveCurlAccountFieldSource())
+ {
+ dSourceVectorFielddExtraDofField += ipv->getConstRefTodSourceVectorFielddExtraDofField()[curlField][extraDOFField];
+ }
+ if (getConstitutiveCurlAccountMecaSource())
+ {
+ dMechFieldSourcedExtraDofField += ipv->getConstRefTodMechanicalFieldSourcedExtraDofField()[curlField][extraDOFField];
+ }
+
+ for (unsigned int j = 0; j < numEdgeShapeFunc; ++j)
+ for (unsigned int k = 0; k < nbFF; ++k)
+ {
+ const unsigned int IndexCurl = 3+getNumNonLocalVariable()+getNumConstitutiveExtraDofDiffusionVariable()+curlField;
+ const unsigned int IndexField = 3+getNumNonLocalVariable()+extraDOFField;
+
+ //mStiff;
+ }
+ }
+ }
+ }
+ }
}
//mStiff.print("bulk");
} // end if sym