Update EffectiveQuantitiesComputer

dune/microstructure/EffectiveQuantitiesComputer.hh

@@ -28,28 +28,28 @@ class EffectiveQuantitiesComputer {
 
 public:
 	static const int dimworld = 3;
-	// static const int nCells = 4;
-    
 	static const int dim = Basis::GridView::dimension;
 
-	using Domain = typename CorrectorComputer<Basis,Material>::Domain; 
 
+	using Domain = typename CorrectorComputer<Basis,Material>::Domain; 
 	using VectorRT = typename CorrectorComputer<Basis,Material>::VectorRT;
 	using MatrixRT = typename CorrectorComputer<Basis,Material>::MatrixRT;
-
 	using Func2Tensor = typename CorrectorComputer<Basis,Material>::Func2Tensor;
 	using FuncVector = typename CorrectorComputer<Basis,Material>::FuncVector;
-
 	using VectorCT = typename CorrectorComputer<Basis,Material>::VectorCT;
-
 	using HierarchicVectorView = typename CorrectorComputer<Basis,Material>::HierarchicVectorView;
+    using Func2int = std::function< int(const Domain&) >;
+    using MatrixPhase = std::function< MatrixRT(const int&) >;
 
 protected:
 
 	CorrectorComputer<Basis,Material>& correctorComputer_; 
-	Func2Tensor prestrain_;
+	// Func2Tensor prestrain_;
+    MatrixPhase prestrain_;
     const Material& material_;
 
+    // Func2int indicatorFunction_;
+
 public:
 	VectorCT B_load_TorusCV_;				//<B, Chi>_L2 
 	// FieldMatrix<double, dim, dim> Q_;	    //effective moduli <LF_i, F_j>_L2
@@ -96,14 +96,21 @@ public:
 	///////////////////////////////
 	// EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer, Func2Tensor prestrain)
     //     : correctorComputer_(correctorComputer), prestrain_(prestrain)
+	// EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer, 
+    //                             Func2Tensor prestrain,
+    //                             const Material& material)
+    //     : correctorComputer_(correctorComputer), 
+    //       prestrain_(prestrain),
+    //       material_(material)
 	EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer, 
-                                Func2Tensor prestrain,
                                 const Material& material)
         : correctorComputer_(correctorComputer), 
-          prestrain_(prestrain),
           material_(material)
     { 
     	
+
+
+        prestrain_ = material_.getPrestrainFunction();
     	// computePrestressLoadCV();
 	  	// computeEffectiveStrains();
         // Q_ = 0;
@@ -156,8 +163,11 @@ public:
                                             correctorComputer_.getCorr_phi3()
 										    };
 
-        auto prestrainGVF  = Dune::Functions::makeGridViewFunction(prestrain_, basis.gridView());
-        auto prestrainFunctional = localFunction(prestrainGVF);   
+        // auto prestrainGVF  = Dune::Functions::makeGridViewFunction(prestrain_, basis.gridView());
+        // auto prestrainFunctional = localFunction(prestrainGVF);   
+
+        auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+        
 
         Q_ = 0 ;
         Bhat_ = 0;
@@ -197,7 +207,8 @@ public:
                 // DerPhi2.bind(e);
                 mu.bind(e);
                 lambda.bind(e);
-                prestrainFunctional.bind(e);
+                // prestrainFunctional.bind(e);
+                localIndicatorFunction.bind(e);
 
                 double elementEnergy = 0.0;
                 double elementPrestrain = 0.0;
@@ -229,12 +240,14 @@ public:
                     auto X1 = G1 + Chi1;
                     //   auto X2 = G2 + Chi2;
                     
-                    
-                    double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, G2);
+                    double energyDensity= scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(G2));
+                    // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, G2);
                     elementEnergy += energyDensity * quadPoint.weight() * integrationElement;      // quad[quadPoint].weight() ???
                     if (b==0)
                     {
-                        elementPrestrain += linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, prestrainFunctional(quadPos)) * quadPoint.weight() * integrationElement;
+                        auto prestrainPhaseValue = prestrain_(localIndicatorFunction(quadPos));
+                        elementPrestrain += scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(prestrainPhaseValue));
+                        // elementPrestrain += linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, prestrainFunctional(quadPos)) * quadPoint.weight() * integrationElement;
                     }
                 }
                 energy += elementEnergy;

src/Cell-Problem-New.cc

@@ -591,7 +591,8 @@ int main(int argc, char *argv[])
 
  
     //--- compute effective quantities
-    auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term,material_);
+    // auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term,material_);
+    auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material_);
     effectiveQuantitiesComputer.computeEffectiveQuantities();