32 #include "young_laplace.h"
35 #include "meshes/simple_rectangular_quadmesh.h"
39 using namespace oomph;
73 Vector<double>& spine_B,
74 Vector< Vector<double> >& dspine_B)
81 dspine_B[0][0] = 1.0 ;
82 dspine_B[1][0] = 0.0 ;
83 dspine_B[0][1] = 0.0 ;
84 dspine_B[1][1] = 1.0 ;
85 dspine_B[0][2] = 0.0 ;
86 dspine_B[1][2] = 0.0 ;
96 double Alpha_min = MathematicalConstants::Pi/2.0*1.5;
99 double Alpha_max = MathematicalConstants::Pi/2.0*0.5;
105 Vector<double>& spine,
106 Vector< Vector<double> >& dspine)
138 template<
class ELEMENT>
154 double new_kappa=
Kappa_pt->value(0)-0.5;
163 void doc_solution(DocInfo& doc_info, ofstream& trace_file);
179 template<
class ELEMENT>
199 Problem::mesh_pt()=
new SimpleRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);
204 if ((n_x%2!=0)||(n_y%2!=0))
206 cout <<
"n_x n_y should be even" << endl;
211 ELEMENT* prescribed_height_element_pt=
dynamic_cast<ELEMENT*
>(
212 mesh_pt()->element_pt(n_y*n_x/2+n_x/2));
215 Control_node_pt=
static_cast<Node*
>(prescribed_height_element_pt->node_pt(0));
217 std::cout <<
"Controlling height at (x,y) : (" << Control_node_pt->x(0)
218 <<
"," << Control_node_pt->x(1) <<
")" <<
"\n" << endl;
222 HeightControlElement* height_control_element_pt=
new HeightControlElement(
226 Kappa_pt=height_control_element_pt->kappa_pt();
241 unsigned n_bound = mesh_pt()->nboundary();
242 for(
unsigned b=0;b<n_bound;b++)
248 unsigned n_node = mesh_pt()->nboundary_node(b);
249 for (
unsigned n=0;n<n_node;n++)
251 mesh_pt()->boundary_node_pt(b,n)->pin(0);
261 unsigned nelement = mesh_pt()->nelement();
262 for(
unsigned i=0;i<nelement;i++)
265 ELEMENT *el_pt =
dynamic_cast<ELEMENT*
>(mesh_pt()->element_pt(i));
277 mesh_pt()->add_element_pt(height_control_element_pt);
280 cout <<
"\nNumber of equations: " << assign_eqn_numbers() << endl;
290 template<
class ELEMENT>
292 ofstream& trace_file)
297 trace_file << -1.0*
Kappa_pt->value(0) <<
" ";
299 trace_file << Control_node_pt->value(0) ;
309 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
311 some_file.open(filename);
312 mesh_pt()->output(some_file,npts);
328 doc_info.set_directory(
"RESLT");
333 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
334 trace_file.open(filename);
338 <<
"VARIABLES=\"<GREEK>k</GREEK>\",\"<GREEK>k</GREEK>_{ex}\",\"h\""
340 trace_file <<
"ZONE" << std::endl;
358 double increment=0.1;
362 for (
unsigned istep=0;istep<nstep;istep++)
369 problem.newton_solve();
2D YoungLaplace problem on rectangular domain, discretised with 2D QYoungLaplace elements....
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc the solution. DocInfo object stores flags/labels for where the output gets written to and the tra...
Node * Control_node_pt
Node at which the height (displacement along spine) is controlled/doced.
YoungLaplaceProblem()
Constructor:
Data * Kappa_pt
Pointer to Data object that stores the prescribed curvature.
void actions_after_newton_solve()
Update the problem after solve: Empty.
void actions_before_newton_solve()
Update the problem before solve.
~YoungLaplaceProblem()
Destructor (empty)
Namespace for "global" problem parameters.
double Alpha_max
Max. spine angle against horizontal plane.
double Controlled_height
Height control value.
double get_exact_kappa()
Exact kappa.
void spine_function(const Vector< double > &x, Vector< double > &spine, Vector< Vector< double > > &dspine)
Spine: The spine vector field as a function of the two coordinates x_1 and x_2, and its derivatives w...
Data * Kappa_pt
Pointer to Data object that stores the prescribed curvature.
void spine_base_function(const Vector< double > &x, Vector< double > &spine_B, Vector< Vector< double > > &dspine_B)
Spine basis: The position vector to the basis of the spine as a function of the two coordinates x_1 a...
double Alpha_min
Min. spine angle against horizontal plane.