From: frohne Date: Sat, 9 Feb 2013 20:46:27 +0000 (+0000) Subject: implement a new class Input to read in data for an obstacle like the Chinese symbol... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=16c2ad6f1c731eaf817e1a97060ba64ccdf78241;p=dealii-svn.git implement a new class Input to read in data for an obstacle like the Chinese symbol for force git-svn-id: https://svn.dealii.org/trunk@28296 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 76a9f13ead..15adf63b26 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -70,6 +70,8 @@ namespace Step42 // @sect3{The PlasticityContactProblem class template} + template class Input; + // This class provides an interface // for a constitutive law. In this // example we are using an elastic @@ -78,8 +80,6 @@ namespace Step42 template class ConstitutiveLaw; - // @sect3{The PlasticityContactProblem class template} - // This class supplies all function // and variables needed to describe // the nonlinear contact problem. It is @@ -99,7 +99,7 @@ namespace Step42 class PlasticityContactProblem { public: - PlasticityContactProblem (int _n_refinements_global, int _n_refinements_local); + PlasticityContactProblem (int _n_refinements_global); void run (); private: @@ -117,7 +117,6 @@ namespace Step42 void output_results (const std::string &title) const; int n_refinements_global; - int n_refinements_local; MPI_Comm mpi_communicator; @@ -151,7 +150,8 @@ namespace Step42 TrilinosWrappers::PreconditionAMG preconditioner_u; TrilinosWrappers::PreconditionAMG preconditioner_t; - std::auto_ptr > plast_lin_hard; + std::unique_ptr > input_obstacle; + std::unique_ptr > plast_lin_hard; double sigma_0; // Yield stress double gamma; // Parameter for the linear isotropic hardening @@ -159,6 +159,199 @@ namespace Step42 double nu; // Poisson ratio }; + template + class Input + { + public: + Input (const char* _name) : + name (_name), + mpi_communicator (MPI_COMM_WORLD), + pcout (std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), + HV (NULL), + lx (0), + ly (0), + nx (0), + ny (0) + {read_surface (name);} + + double hv(int i, int j); + + double& set_height(int i, int j); + + double mikro_height(double x,double y, double z); + + void read_surface(const char* name); + + private: + const char* name; + MPI_Comm mpi_communicator; + ConditionalOStream pcout; + double* HV; + double lx, ly; + double hx, hy; + int nx, ny; + }; + + template + double Input::hv(int i, int j) { + return HV[nx*j+i]; // i indiziert x-werte, j indiziert y-werte + } + + template + double& Input::set_height(int i, int j) { + return HV[nx*j+i]; // i indiziert x-werte, j indiziert y-werte + } + + template + double Input::mikro_height(double x,double y, double z) { + int ix = (int)(x/hx); + int iy = (int)(y/hy); + + if (ix<0) { + ix = 0; + // cerr << "hm\n"; + } + if (iy<0) { + iy = 0; + // cerr << "hm\n"; + } + + if (ix>=nx-1) { + ix = nx-2; + } + if (iy>=ny-1) { + iy = ny-2; + } + + double val = 0.; + { + FullMatrix H(4,4); + Vector X(4); + Vector b(4); + + double xx = 0.; + double yy = 0.; + + xx = (ix )*hx; + yy = (iy )*hy; + H(0,0) = xx; + H(0,1) = yy; + H(0,2) = xx*yy; + H(0,3) = 1.; + b(0) = hv(ix ,iy ); + + xx = (ix+1)*hx; + yy = (iy )*hy; + H(1,0) = xx; + H(1,1) = yy; + H(1,2) = xx*yy; + H(1,3) = 1.; + b(1) = hv(ix+1,iy ); + + xx = (ix+1)*hx; + yy = (iy+1)*hy; + H(2,0) = xx; + H(2,1) = yy; + H(2,2) = xx*yy; + H(2,3) = 1.; + b(2) = hv(ix+1,iy+1); + + xx = (ix )*hx; + yy = (iy+1)*hy; + H(3,0) = xx; + H(3,1) = yy; + H(3,2) = xx*yy; + H(3,3) = 1.; + b(3) = hv(ix ,iy+1); + + H.gauss_jordan(); + H.vmult(X,b); + + val = X(0)*x + X(1)*y + X(2)*x*y + X(3); + } + + return val; + } + + template + void Input::read_surface(const char* name) { + int SZ = 100000; + FILE* fp = fopen(name,"r"); + char* zeile = new char[SZ]; + char* hlp_str = new char[SZ]; + + double hlp; + + int POS; + //////////////////////////////// + fgets(zeile,SZ,fp); + POS = strcspn(zeile,"="); + for (int i=0;i<=POS;i++) { + zeile[i] = ' '; + } + sscanf(zeile,"%d",&nx); + //////////////////////////////// + fgets(zeile,SZ,fp); + POS = strcspn(zeile,"="); + for (int i=0;i<=POS;i++) { + zeile[i] = ' '; + } + sscanf(zeile,"%d",&ny); + //////////////////////////////// + fgets(zeile,SZ,fp); + POS = strcspn(zeile,"="); + for (int i=0;i<=POS;i++) { + zeile[i] = ' '; + } + sscanf(zeile,"%lf",&lx); + //////////////////////////////// + fgets(zeile,SZ,fp); + POS = strcspn(zeile,"="); + for (int i=0;i<=POS;i++) { + zeile[i] = ' '; + } + sscanf(zeile,"%lf",&ly); + + pcout<< nx << " " << ny << " " << lx << " " << ly << " " < max_hlp) + max_hlp=hlp; + if (hlp < min_hlp) + min_hlp=hlp; + + set_height(k,ny-1-j) = hlp; + int pos = strcspn(zeile,","); + if (!strpbrk(zeile,",")) { + reached = 1; + continue; + } + for (int i=0;i<=pos;i++) { + zeile[i] = ' '; + } + } + j++; + } + pcout<< "/** highest point: " << max_hlp < class ConstitutiveLaw { @@ -379,13 +572,18 @@ namespace Step42 class Obstacle : public Function { public: - Obstacle () : Function(dim) {}; + Obstacle (std::unique_ptr > const &_input) : + Function(dim), + input_obstacle_copy(std::move (_input)) {}; virtual double value (const Point &p, const unsigned int component = 0) const; virtual void vector_value (const Point &p, Vector &values) const; + + private: + std::unique_ptr > const &input_obstacle_copy; }; template @@ -393,7 +591,7 @@ namespace Step42 const unsigned int component) const { double R = 0.03; - double return_value = 0.0; + double return_value = 100.0; if (component == 0) return_value = p(0); if (component == 1) @@ -401,12 +599,19 @@ namespace Step42 if (component == 2) { // Hindernis Dortmund - double x1 = p(0); - double x2 = p(1); - if (((x2-0.5)*(x2-0.5)+(x1-0.5)*(x1-0.5)<=0.3*0.3)&&((x2-0.5)*(x2-0.5)+(x1-1.0)*(x1-1.0)>=0.4*0.4)&&((x2-0.5)*(x2-0.5)+x1*x1>=0.4*0.4)) - return_value = 0.999; - else - return_value = 1e+10; +// double x1 = p(0); +// double x2 = p(1); +// if (((x2-0.5)*(x2-0.5)+(x1-0.5)*(x1-0.5)<=0.3*0.3)&&((x2-0.5)*(x2-0.5)+(x1-1.0)*(x1-1.0)>=0.4*0.4)&&((x2-0.5)*(x2-0.5)+x1*x1>=0.4*0.4)) +// return_value = 0.999; +// else +// return_value = 1e+10; + + // Hindernis Werkzeug TKSE + return_value = 1.999 - input_obstacle_copy->mikro_height (p(0), p(1), p(2)); +// std::cout<< "Obstacle value: " << return_value +// << " p(0) = " << p(0) +// << " p(1) = " << p(1) +// < - PlasticityContactProblem::PlasticityContactProblem (int _n_refinements_global, int _n_refinements_local) + PlasticityContactProblem::PlasticityContactProblem (int _n_refinements_global) : n_refinements_global (_n_refinements_global), - n_refinements_local (_n_refinements_local), mpi_communicator (MPI_COMM_WORLD), triangulation (mpi_communicator), fe (FE_Q(1), dim), @@ -860,7 +1064,7 @@ namespace Step42 { clock_t start_proj, end_proj; - const EquationData::Obstacle obstacle; + const EquationData::Obstacle obstacle (input_obstacle); std::vector vertex_touched (dof_handler.n_dofs (), false); typename DoFHandler::active_cell_iterator @@ -1032,7 +1236,7 @@ namespace Step42 PrimitiveVectorMemory mem; TrilinosWrappers::MPI::Vector tmp (system_rhs_newton); - const double solver_tolerance = 1e-4 * + const double solver_tolerance = 1e-3 * system_matrix_newton.residual (tmp, distributed_solution, system_rhs_newton); SolverControl solver_control (system_matrix_newton.m(), solver_tolerance); @@ -1314,7 +1518,10 @@ namespace Step42 Timer t; run_time.resize (8); - const unsigned int n_cycles = 5; + // Read in the obstacle data. + input_obstacle.reset (new Input("obstacle_file.dat")); + + const unsigned int n_cycles = 6; for (unsigned int cycle=0; cycle laplace_problem_3d (_n_refinements_global, _n_refinements_local); + PlasticityContactProblem<3> laplace_problem_3d (_n_refinements_global); laplace_problem_3d.run (); }