diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc index 7fa988ff6d90f309a6df913a84dde22b49d02d8f..38f9d535e7ca4d8a1f7ac4fac431664f5a06a085 100644 --- a/src/dune-microstructure.cc +++ b/src/dune-microstructure.cc @@ -551,7 +551,7 @@ void computeElementLoadVector( const LocalView& localView, double energyDensity = 0.0; for (int ii=0; ii<nCompo; ii++) for (int jj=0; jj<nCompo; jj++) - energyDensity += stressV[ii][jj] * forceTerm(quadPos)[ii][jj]; + energyDensity += stressV[ii][jj] *forceTerm(quadPos)[ii][jj]; //TEST VERWENDE NEGATIVE DER LAST size_t row = localView.tree().child(k).localIndex(i); @@ -577,7 +577,7 @@ void computeElementLoadVector( const LocalView& localView, double energyDensityfG = 0; for (int ii=0; ii<nCompo; ii++) for (int jj=0; jj<nCompo; jj++) - energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj]; + energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj]; //TEST VERWENDE NEGATIVE DER LAST elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement; } @@ -1635,8 +1635,7 @@ int main(int argc, char *argv[]) /////////////////////////////////// // Get Parameters/Data /////////////////////////////////// - - double gamma = parameterSet.get<double>("gamma",2.0); // ratio dimension reduction to homogenization + double gamma = parameterSet.get<double>("gamma",2000.0); // ratio dimension reduction to homogenization // Material parameter laminat // double E1 = parameterSet.get<double>("E1", 17e6); //material1 @@ -1645,8 +1644,6 @@ int main(int argc, char *argv[]) // double E2 = parameterSet.get<double>("E2", 35e6); //material2 // double nu2 = parameterSet.get<double>("nu2", 0.5); - - // // Plate geometry settings double width = parameterSet.get<double>("width", 1.0); //geometry cell, cross section // double len = parameterSet.get<double>("len", 10.0); //length @@ -1676,9 +1673,6 @@ int main(int argc, char *argv[]) log << "alpha: " << alpha << std::endl; log << "gamma: " << gamma << std::endl; - - - /////////////////////////////////// // --- Choose Solver --- @@ -1694,7 +1688,7 @@ int main(int argc, char *argv[]) /////////////////////////////////// // Generate the grid /////////////////////////////////// - using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim> >; + using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >; // Domain : [0,1)^2 x (-1/2, 1/2) // FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0}); @@ -1704,25 +1698,22 @@ int main(int argc, char *argv[]) FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0}); FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0}); - int nE = 10; + int nE = 20; log << "Number of Elements: " << nE << std::endl; debugLog << "Number of Elements: " << nE << std::endl; std::array<int,dim> nElements={nE,nE,nE}; +// std::array<int,dim> nElements={100,100,1}; // TEST // printvector(std::cout, coeffT_1 , "coeffT_1", "--" ); // printvector(std::cout, x_1 , "x_1", "--" );; //#Elements in each direction - CellGridType grid_CE(lower,upper,nElements); //Corrector problem Domain: + CellGridType grid_CE(lower,upper,nElements); //Corrector Problem Domain: using GridView = CellGridType::LeafGridView; const GridView gridView_CE = grid_CE.leafGridView(); - - - - using ScalarRT = FieldVector< double, 1>; using VectorRT = FieldVector< double, nCompo>; using MatrixRT = FieldMatrix< double, nCompo, nCompo>; @@ -1730,7 +1721,6 @@ int main(int argc, char *argv[]) using FuncScalar = std::function< ScalarRT(const Domain&) >; using Func2Tensor = std::function< MatrixRT(const Domain&) >; - using VectorCT = BlockVector<FieldVector<double,1> >; using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >; // using VectorCT = BlockVector<FieldVector<double,dim> >; @@ -1770,8 +1760,6 @@ int main(int argc, char *argv[]) else return mu2; // return mu1; - - }; @@ -1783,9 +1771,6 @@ int main(int argc, char *argv[]) }; - - - auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE); auto muLocal = localFunction(muGridF); auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE); @@ -1798,8 +1783,6 @@ int main(int argc, char *argv[]) - - ///////////////////////////////////////////////////////// // Choose a finite element space for Cell Problem ///////////////////////////////////////////////////////// @@ -1948,6 +1931,19 @@ int main(int argc, char *argv[]) assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1); assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2); assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3); + + + // TEST multiply RHS with -1 : + for(int i=0; i< Basis_CE.size()+3 ; i++) + { + load_alpha1[i] = (-1.0)*load_alpha1[i]; + load_alpha2[i] = (-1.0)*load_alpha2[i]; + load_alpha3[i] = (-1.0)*load_alpha3[i]; + } + + + + // printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");