oomph-lib

00001 //LIC// ====================================================================
00002 //LIC// This file forms part of oomph-lib, the object-oriented, 
00003 //LIC// multi-physics finite-element library, available 
00004 //LIC// at http://www.oomph-lib.org.
00005 //LIC// 
00006 //LIC//           Version 0.90. August 3, 2009.
00007 //LIC// 
00008 //LIC// Copyright (C) 2006-2009 Matthias Heil and Andrew Hazel
00009 //LIC// 
00010 //LIC// This library is free software; you can redistribute it and/or
00011 //LIC// modify it under the terms of the GNU Lesser General Public
00012 //LIC// License as published by the Free Software Foundation; either
00013 //LIC// version 2.1 of the License, or (at your option) any later version.
00014 //LIC// 
00015 //LIC// This library is distributed in the hope that it will be useful,
00016 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 //LIC// Lesser General Public License for more details.
00019 //LIC// 
00020 //LIC// You should have received a copy of the GNU Lesser General Public
00021 //LIC// License along with this library; if not, write to the Free Software
00022 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
00023 //LIC// 02110-1301  USA.
00024 //LIC// 
00025 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
00026 //LIC// 
00027 //LIC//====================================================================
00028 // Demonstrate use of geometric object as GeneralisedElement
00029 
00030  
00031 // Generic oomph-lib headers
00032 #include "generic.h"
00033 
00034 // Circle as generalised element:
00035 #include "circle_as_generalised_element.h"
00036 
00037 using namespace std;
00038 
00039 using namespace oomph;
00040 
00041 ///////////////////////////////////////////////////////////////////////
00042 ///////////////////////////////////////////////////////////////////////
00043 ///////////////////////////////////////////////////////////////////////
00044 
00045 
00046 //======start_of_problem==============================================
00047 /// Problem to demonstrate the use of a GeomObject as a 
00048 /// GeneralisedElement: A geometric object (a Circle) is "upgraded"
00049 /// to a GeneralisedElement. The position of the Circle is
00050 /// determined by a balance of forces, assuming that the
00051 /// Circle is mounted on an elastic spring of specified
00052 /// stiffness and loaded by a vertical "load". 
00053 //====================================================================
00054 class GeomObjectAsGeneralisedElementProblem : public Problem
00055 {
00056 
00057 public:
00058 
00059  /// Constructor
00060  GeomObjectAsGeneralisedElementProblem();
00061 
00062  /// Update the problem specs after solve (empty)
00063  void actions_after_newton_solve(){}
00064   
00065  /// Update the problem specs before solve (empty)
00066  void actions_before_newton_solve() {}
00067 
00068  /// Doc the solution
00069  void doc_solution();
00070 
00071  /// Return value of the "load" on the elastically supported ring
00072  double& load()
00073   {
00074    return *Load_pt->value_pt(0);
00075   }
00076 
00077  /// Access to DocInfo object
00078  DocInfo& doc_info() {return Doc_info;}
00079 
00080 private:
00081 
00082  /// Trace file
00083  ofstream Trace_file;
00084 
00085  /// \short Pointer to data item that stores the "load" on the ring
00086  Data* Load_pt;
00087 
00088  /// Doc info object
00089  DocInfo Doc_info;
00090 
00091 };
00092 
00093 
00094 
00095 
00096 
00097 //=============================start_of_problem_constructor===============
00098 /// Constructor
00099 //========================================================================
00100 GeomObjectAsGeneralisedElementProblem::GeomObjectAsGeneralisedElementProblem()
00101 { 
00102  
00103  // Set coordinates and radius for the circle
00104  double x_c=0.5;
00105  double y_c=0.0;
00106  double R=1.0;
00107 
00108  // Build GeomObject that's been upgraded to a GeneralisedElement
00109  // GeneralisedElement* 
00110  ElasticallySupportedRingElement* geom_object_element_pt = 
00111   new ElasticallySupportedRingElement(x_c,y_c,R);
00112 
00113  // Set the stiffness of the elastic support
00114  geom_object_element_pt->k_stiff()=0.3;
00115 
00116  // Build mesh
00117  mesh_pt()=new Mesh;
00118 
00119  // So far, the mesh is completely empty. Let's add the 
00120  // one (and only!) GeneralisedElement to it:
00121  mesh_pt()->add_element_pt(geom_object_element_pt);
00122 
00123  // Create the load (a Data object with a single value)
00124  Load_pt=new Data(1);
00125    
00126  // The load is prescribed so its one-and-only value is pinned
00127  Load_pt->pin(0);
00128 
00129  // Set the pointer to the Data object that specifies the 
00130  // load on the ring
00131  geom_object_element_pt->set_load_pt(Load_pt);
00132  
00133  // Setup equation numbering scheme.
00134  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl; 
00135 
00136  // Set output directory
00137  Doc_info.set_directory("RESLT"); 
00138   
00139  // Open trace file
00140  char filename[100];
00141  sprintf(filename,"%s/trace.dat",Doc_info.directory().c_str());
00142  Trace_file.open(filename);
00143  Trace_file << "VARIABLES=\"load\",\"y<sub>circle</sub>\"" << std::endl;
00144 
00145 } // end of constructor
00146 
00147 
00148 
00149 
00150 //===========================start_of_doc_solution========================
00151 /// Doc the solution in tecplot format.
00152 //========================================================================
00153 void GeomObjectAsGeneralisedElementProblem::doc_solution()
00154 { 
00155 
00156  ofstream some_file;
00157  char filename[100];
00158 
00159  // Number of plot points
00160  unsigned npts=100;
00161 
00162  // Lagrangian coordinate and position vector (both as vectors)
00163  Vector<double> zeta(1);
00164  Vector<double> r(2);
00165  
00166  // Output solution 
00167  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
00168          Doc_info.number());
00169  some_file.open(filename);
00170  for (unsigned i=0;i<npts;i++)
00171   {
00172    zeta[0]=2.0*MathematicalConstants::Pi*double(i)/double(npts-1);
00173    static_cast<ElasticallySupportedRingElement*>(mesh_pt()->element_pt(0))->
00174     position(zeta,r);  
00175    some_file << r[0] << " " << r[1] << std::endl;
00176   }
00177  some_file.close();
00178 
00179 
00180  // Write "load" and vertical position of the ring's centre
00181  Trace_file 
00182   << static_cast<ElasticallySupportedRingElement*>(
00183    mesh_pt()->element_pt(0))->load()
00184   << " "
00185   << static_cast<ElasticallySupportedRingElement*>(
00186    mesh_pt()->element_pt(0))->y_c()
00187   << " "
00188   << std::endl;
00189 
00190 } // end of doc_solution
00191 
00192 
00193  
00194 
00195 
00196 
00197 
00198 
00199 ////////////////////////////////////////////////////////////////////////
00200 ////////////////////////////////////////////////////////////////////////
00201 ////////////////////////////////////////////////////////////////////////
00202 
00203 
00204 
00205 
00206 
00207 //===============start_of_driver==========================================
00208 /// Driver 
00209 //========================================================================
00210 int main()
00211 {
00212 
00213  // Set up the problem 
00214  GeomObjectAsGeneralisedElementProblem problem;
00215  
00216  // Initial value for the load 
00217  problem.load()=-0.3;
00218  
00219  // Loop for different loads
00220  //-------------------------
00221  
00222  // Number of steps
00223  unsigned nstep=5;
00224  
00225  // Increment in load
00226  double dp=0.6/double(nstep-1);
00227  
00228  for (unsigned istep=0;istep<nstep;istep++)
00229   {    
00230    // Solve/doc
00231    problem.newton_solve();
00232    problem.doc_solution();
00233    
00234    //Increment counter for solutions 
00235    problem.doc_info().number()++;    
00236    
00237    // Change load on ring
00238    problem.load()+=dp;
00239   } 
00240  
00241 } // end of driver
00242 
00243
geom_object_element.cc
