direct sparse solver & umicore compatibility (Martin Horak)

eudoxos · Apr 25, 2017 · 5b9bbfa · 5b9bbfa
1 parent eeb051f
commit 5b9bbfa
Show file tree

Hide file tree

Showing 14 changed files with 62 additions and 23 deletions.
diff --git a/bindings/python/oofemlib.cpp b/bindings/python/oofemlib.cpp
@@ -699,7 +699,7 @@ void pyclass_Element()
         .add_property("length", &PyElement::computeLength)
         .def("giveLabel", &Element::giveLabel)
         .add_property("label",&Element::giveLabel)
-        .def("giveLocationArray", &PyElement::giveLocationArray)
+        //.def("giveLocationArray", &PyElement::giveLocationArray)
         .def("giveNumberOfDofManagers", &PyElement::giveNumberOfDofManagers)
         .add_property("numberOfDofManagers", &PyElement::giveNumberOfDofManagers)
         .def("computeNumberOfDofs", &PyElement::computeNumberOfDofs)
@@ -810,7 +810,7 @@ void pyclass_Load()
 {
     class_<PyLoad, bases<GeneralBoundaryCondition>, boost::noncopyable >("Load", no_init)
         .def("setComponentArray", &Load::setComponentArray)
-        .def("giveCopyOfComponentArray", &Load::giveCopyOfComponentArray)
+        .def("giveComponentArray", &Load::giveCopyOfComponentArray, return_value_policy<manage_new_object>())
         .def("computeValueAt", pure_virtual( &Load::computeValueAt))
         ;
 }

diff --git a/src/dss/SparseGridMtxLL.cpp b/src/dss/SparseGridMtxLL.cpp
@@ -367,7 +367,7 @@ void SparseGridMtxLL :: Factorize()
             BlockArith->LL_Decomposition(dd + Djj);
 
             if ( ( bj % newdotcnount ) == 0 ) {
-                Write(".");
+                //Write(".");
             }
 
             Djj += block_storage;
@@ -471,7 +471,7 @@ void SparseGridMtxLL :: Factorize_Incomplete()
             BlockArith->LL_Decomposition(dd + Djj);
 
             if ( ( bj % newdotcnount ) == 0 ) {
-                Write(".");
+                //Write(".");
             }
 
             Djj += block_storage;

diff --git a/src/dss/SparseGridMtxLU.cpp b/src/dss/SparseGridMtxLU.cpp
@@ -377,7 +377,7 @@ void SparseGridMtxLU :: Factorize()
             BlockArith->LU_Decomposition(dd + Djj);
 
             if ( ( bj % newdotcnount ) == 0 ) {
-                Write(".");
+                //Write(".");
             }
 
             Djj += block_storage;

diff --git a/src/oofemlib/cltypes.C b/src/oofemlib/cltypes.C
@@ -96,6 +96,7 @@ InternalStateValueType giveInternalStateValueType(InternalStateType type)
     case IST_DeviatoricStrain:
     case IST_DeviatoricStress:
     case IST_CauchyStressTensor:
+    case IST_StrainTensorStressDependent:
         return ISVT_TENSOR_S3;
 
     case IST_BeamForceMomentumTensor:
@@ -126,6 +127,7 @@ InternalStateValueType giveInternalStateValueType(InternalStateType type)
     case IST_CrackDirs:
     case IST_CrackStatuses:
     case IST_CrackVector:
+    case IST_PrincipalStrainTensorStressDependent:
         return ISVT_VECTOR;
 
     case IST_MaxEquivalentStrainLevel:

diff --git a/src/oofemlib/internalstatetype.h b/src/oofemlib/internalstatetype.h
@@ -149,7 +149,9 @@ namespace oofem {
     ENUM_ITEM_WITH_VALUE(IST_InterfaceFirstPKTraction, 99) \
     ENUM_ITEM_WITH_VALUE(IST_StressTensor_Reduced, 100) \
     ENUM_ITEM_WITH_VALUE(IST_StrainTensor_Reduced, 101) \
-    ENUM_ITEM_WITH_VALUE(IST_CrossSectionNumber, 102)
+    ENUM_ITEM_WITH_VALUE(IST_CrossSectionNumber, 102) \
+    ENUM_ITEM_WITH_VALUE(IST_StrainTensorStressDependent, 103) \
+    ENUM_ITEM_WITH_VALUE(IST_PrincipalStrainTensorStressDependent, 104)
 /**
  * Type  representing the physical meaning of element or constitutive model internal variable.
  * Values of this type are used, when these internal variables are requested.

diff --git a/src/oofemlib/load.h b/src/oofemlib/load.h
@@ -165,8 +165,7 @@ class OOFEM_EXPORT Load : public GeneralBoundaryCondition
 
 public:
     void setComponentArray(FloatArray &arry) { componentArray = arry; }
-    FloatArray giveCopyOfComponentArray() { FloatArray answer = componentArray;
-                                            return answer; }
+    FloatArray* giveCopyOfComponentArray() { return new FloatArray(componentArray); }
 };
 } // end namespace oofem
 #endif // load_h
diff --git a/src/oofemlib/nrsolver.C b/src/oofemlib/nrsolver.C
@@ -275,6 +275,7 @@ NRSolver :: solve(SparseMtrx *k, FloatArray *R, FloatArray *R0,
             break;
         } else if ( nite >= nsmax ) {
             OOFEM_LOG_DEBUG("Maximum number of iterations reached\n");
+				OOFEM_ERROR("ASG\n");
             break;
         }
 

diff --git a/src/oofemlib/zznodalrecoverymodel.C b/src/oofemlib/zznodalrecoverymodel.C
@@ -133,7 +133,7 @@ ZZNodalRecoveryModel :: recoverValues(InternalStateType type, TimeStep *tStep)
                 rhs.resize(regionDofMans, regionValSize);
                 rhs.zero();
                 if ( regionValSize == 0 ) {
-                    OOFEM_LOG_RELEVANT( "ZZNodalRecoveryModel :: unknown size of InternalStateType %s\n", __InternalStateTypeToString(type) );
+                    //OOFEM_LOG_RELEVANT( "ZZNodalRecoveryModel :: unknown size of InternalStateType %s\n", __InternalStateTypeToString(type) );
                 }
             } else if ( regionValSize != nsig.giveNumberOfColumns() ) {
                 nsig.resize(regionDofMans, regionValSize);

diff --git a/src/sm/interfaceelement1d.C b/src/sm/interfaceelement1d.C
@@ -246,6 +246,7 @@ InterfaceElem1d :: computeVolumeAround(GaussPoint *gp)
 // Returns the length of the receiver. This method is valid only if 1
 // Gauss point is used.
 {
+	 return this->giveCrossSection()->give(CS_Area,gp);
     return 1.0; //MUST be set to 1.0
 }
 

diff --git a/src/sm/nldeidynamic.C b/src/sm/nldeidynamic.C
@@ -78,7 +78,7 @@ NumericalMethod *NlDEIDynamic :: giveNumericalMethod(MetaStep *mStep)
 //     - SolutionOfLinearEquations
 
 {
-    //return NULL;  // Not necessary here - Diagonal matrix and simple inversion is used.
+    return NULL;  // Not necessary here - Diagonal matrix and simple inversion is used.
     if ( nMethod ) {
         return nMethod;
     }
@@ -618,7 +618,7 @@ NlDEIDynamic :: computeMassMtrx(FloatArray &massMatrix, double &maxOm, TimeStep
         for ( j = 1; j <= n; j++ ) {
             jj = loc.at(j);
             if ( ( jj ) && ( charMtrx.at(j, j) <= maxElmass * ZERO_REL_MASS ) ) {
-                charMtrx.at(j, j) = charMtrx2.at(j, j) / maxOmEl;
+					if (maxOmEl > 0.) charMtrx.at(j, j) = charMtrx2.at(j, j) / maxOmEl;
             }
         }
 #endif

diff --git a/src/sm/rankinemat.C b/src/sm/rankinemat.C
@@ -169,6 +169,8 @@ RankineMat :: CreateStatus(GaussPoint *gp) const
 void
 RankineMat :: giveRealStressVector_1d(FloatArray &answer, GaussPoint *gp, const FloatArray &totalStrain, TimeStep *tStep)
 {
+	 FloatArray reducedTotalStrainVector;
+    this->giveStressDependentPartOfStrainVector(reducedTotalStrainVector, gp, totalStrain, tStep, VM_Total);
     FloatArray strainVector;
     FloatMatrix d;
     RankineMatStatus *status = static_cast< RankineMatStatus * >( this->giveStatus(gp) );
@@ -178,7 +180,7 @@ RankineMat :: giveRealStressVector_1d(FloatArray &answer, GaussPoint *gp, const
     this->initGpForNewStep(gp);
 
     // elastoplasticity
-    this->performPlasticityReturn(gp, totalStrain);
+    this->performPlasticityReturn(gp, reducedTotalStrainVector);
 
     // damage
     double omega = computeDamage(gp, tStep);
@@ -202,14 +204,16 @@ RankineMat :: giveRealStressVector_PlaneStress(FloatArray &answer,
                                                const FloatArray &totalStrain,
                                                TimeStep *tStep)
 {
+	 FloatArray reducedTotalStrainVector;
+    this->giveStressDependentPartOfStrainVector(reducedTotalStrainVector, gp, totalStrain, tStep, VM_Total);
     RankineMatStatus *status = static_cast< RankineMatStatus * >( this->giveStatus(gp) );
 
     // initialization
     this->initTempStatus(gp);
     this->initGpForNewStep(gp);
 
     // elastoplasticity
-    this->performPlasticityReturn(gp, totalStrain);
+    this->performPlasticityReturn(gp, reducedTotalStrainVector);
 
     // damage
     double omega = computeDamage(gp, tStep);

diff --git a/src/sm/simpleinterfacemat.C b/src/sm/simpleinterfacemat.C
@@ -122,14 +122,14 @@ SimpleInterfaceMaterial :: giveRealStressVector(FloatArray &answer, GaussPoint *
     case _2dInterface:
         answer.resize(2);
         shearStrain.at(1) = strainVector.at(2);
-        shearStress.beScaled(this->kn, shearStrain);
+        shearStress.beScaled(this->ks, shearStrain);
         shearStress.subtract(tempShearStressShift);
         dp = shearStress.dotProduct(shearStress, 1);
         if ( dp > maxShearStress * maxShearStress ) {
             shearStress.times( maxShearStress / sqrt(dp) );
         }
 
-        tempShearStressShift.beScaled(this->kn, shearStrain);
+        tempShearStressShift.beScaled(this->ks, shearStrain);
         tempShearStressShift.subtract(shearStress);
         answer.at(2) = shearStress.at(1);
         break;
@@ -138,14 +138,14 @@ SimpleInterfaceMaterial :: giveRealStressVector(FloatArray &answer, GaussPoint *
         answer.resize(3);
         shearStrain.at(1) = strainVector.at(2);
         shearStrain.at(2) = strainVector.at(3);
-        shearStress.beScaled(this->kn, shearStrain);
+        shearStress.beScaled(this->ks, shearStrain);
         shearStress.subtract(tempShearStressShift);
         dp = shearStress.dotProduct(shearStress, 2);
         if ( dp > maxShearStress * maxShearStress ) {
             shearStress.times( maxShearStress / sqrt(dp) );
         }
 
-        tempShearStressShift.beScaled(this->kn, shearStrain);
+        tempShearStressShift.beScaled(this->ks, shearStrain);
         tempShearStressShift.subtract(shearStress);
         answer.at(2) = shearStress.at(1);
         answer.at(3) = shearStress.at(2);
@@ -205,13 +205,16 @@ SimpleInterfaceMaterial :: giveStiffnessMatrix(FloatMatrix &answer,
         answer.resize(2, 2);
         if ( rMode == SecantStiffness || rMode == TangentStiffness ) {
             if ( normalStrain + normalClearance <= 0. ) {
-                answer.at(1, 1) = answer.at(2, 2) = this->kn; //in compression and after the clearance gap closed
+                answer.at(1, 1) = this->kn; //in compression and after the clearance gap closed
+                answer.at(2, 2) = this->ks;
             } else {
-                answer.at(1, 1) = answer.at(2, 2) = this->kn * this->stiffCoeff;
+                answer.at(1, 1) = this->kn * this->stiffCoeff;
+                answer.at(2, 2) = this->ks * this->stiffCoeff;
             }
         } else {
             if ( rMode == ElasticStiffness ) {
-                answer.at(1, 1) = answer.at(2, 2) = this->kn;
+                answer.at(1, 1) = this->kn;
+                answer.at(2, 2) = this->ks;
             } else {
                 _error2( "give2dInterfaceMaterialStiffnessMatrix: unknown MatResponseMode (%s)", __MatResponseModeToString(rMode) );
             }
@@ -223,13 +226,16 @@ SimpleInterfaceMaterial :: giveStiffnessMatrix(FloatMatrix &answer,
         answer.resize(3, 3);
         if ( rMode == SecantStiffness || rMode == TangentStiffness ) {
             if ( normalStrain + normalClearance <= 0. ) {
-                answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = this->kn; //in compression and after the clearance gap closed
+                answer.at(1, 1) = this->kn; //in compression and after the clearance gap closed
+                answer.at(2, 2) = answer.at(3, 3) = this->ks;
             } else {
-                answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = this->kn * this->stiffCoeff;
+                answer.at(1, 1) = this->kn * this->stiffCoeff;
+                answer.at(2, 2) = answer.at(3, 3) = this->ks * this->stiffCoeff;
             }
         } else {
             if ( rMode == ElasticStiffness ) {
-                answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = this->kn;
+                answer.at(1, 1) = this->kn;
+                answer.at(2, 2) = answer.at(3, 3) = this->ks;
             } else {
                 _error2( "give2dInterfaceMaterialStiffnessMatrix: unknown MatResponseMode (%s)", __MatResponseModeToString(rMode) );
             }
@@ -274,10 +280,13 @@ SimpleInterfaceMaterial :: initializeFrom(InputRecord *ir)
     frictCoeff = 0.;
     stiffCoeff = 0.;
     normalClearance = 0.;
+	 ks = -1;
     IR_GIVE_FIELD(ir, kn, _IFT_SimpleInterfaceMaterial_kn);
+    IR_GIVE_OPTIONAL_FIELD(ir, ks, _IFT_SimpleInterfaceMaterial_ks);
     IR_GIVE_OPTIONAL_FIELD(ir, frictCoeff, _IFT_SimpleInterfaceMaterial_frictCoeff);
     IR_GIVE_OPTIONAL_FIELD(ir, stiffCoeff, _IFT_SimpleInterfaceMaterial_stiffCoeff);
     IR_GIVE_OPTIONAL_FIELD(ir, normalClearance, _IFT_SimpleInterfaceMaterial_normalClearance);
+	 if (ks<0) ks = ks;
 
     return StructuralMaterial :: initializeFrom(ir);
 }
@@ -288,6 +297,7 @@ SimpleInterfaceMaterial :: giveInputRecord(DynamicInputRecord &input)
 {
     StructuralMaterial :: giveInputRecord(input);
     input.setField(this->kn, _IFT_SimpleInterfaceMaterial_kn);
+    input.setField(this->ks, _IFT_SimpleInterfaceMaterial_ks);
     input.setField(this->frictCoeff, _IFT_SimpleInterfaceMaterial_frictCoeff);
     input.setField(this->stiffCoeff, _IFT_SimpleInterfaceMaterial_stiffCoeff);
     input.setField(this->normalClearance, _IFT_SimpleInterfaceMaterial_normalClearance);

diff --git a/src/sm/simpleinterfacemat.h b/src/sm/simpleinterfacemat.h
@@ -48,6 +48,7 @@
 //@{
 #define _IFT_SimpleInterfaceMaterial_Name "simpleintermat"
 #define _IFT_SimpleInterfaceMaterial_kn "kn"
+#define _IFT_SimpleInterfaceMaterial_ks "ks"
 #define _IFT_SimpleInterfaceMaterial_knt "knt"
 #define _IFT_SimpleInterfaceMaterial_frictCoeff "fc"
 #define _IFT_SimpleInterfaceMaterial_stiffCoeff "stiffcoeff"
@@ -96,6 +97,7 @@ class SimpleInterfaceMaterial : public StructuralMaterial
 {
 protected:
     double kn;
+    double ks;
     double stiffCoeff;
     double frictCoeff;
     /// Normal distance which needs to be closed when interface element should act in compression (distance is 0 by default).

diff --git a/src/sm/structuralmaterial.C b/src/sm/structuralmaterial.C
@@ -1796,6 +1796,24 @@ StructuralMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalSt
     } else if ( type == IST_FirstPKStressTensor ) {
         answer = status->givePVector();
         return 1;
+    } else if ( type == IST_StrainTensorStressDependent ) {
+        ///@todo Fill in correct full form values here! This just adds zeros!
+		  FloatArray s;
+		  this->giveStressDependentPartOfStrainVector(s,gp,status->giveStrainVector(),tStep,VM_Total);
+        StructuralMaterial :: giveFullSymVectorForm( answer, s, gp->giveMaterialMode() );
+        if ( gp->giveMaterialMode() == _PlaneStress ) {
+            double Nxy = this->give(NYxy, gp);
+            double Nxz = this->give(NYxz, gp);
+            double Nyz = this->give(NYyz, gp);
+            double Nyx = Nxy * this->give(Ey, gp) / this->give(Ex, gp);
+            answer.at(3) = ( -( Nxz + Nxy * Nyz ) * answer.at(1) - ( Nyz + Nxz * Nyx ) * answer.at(2) ) / ( 1. - Nxy * Nyx );
+        }
+        return 1;
+    } else if ( type == IST_PrincipalStrainTensorStressDependent ) {
+        FloatArray s;
+		  this->giveStressDependentPartOfStrainVector(s,gp,status->giveStrainVector(),tStep,VM_Total);
+        this->computePrincipalValues(answer, s, principal_strain);
+        return 1;
     } else {
         return Material :: giveIPValue(answer, gp, type, tStep);
     }
-Original file line number
+Diff line change
@@ Expand Up @@
                 break;
             } else if ( nite >= nsmax ) {
                 OOFEM_LOG_DEBUG("Maximum number of iterations reached\n");
+    				OOFEM_ERROR("ASG\n");
                 break;
             }
@@ Expand Down @@