32 #ifndef OOMPH_PROJECTION_HEADER 33 #define OOMPH_PROJECTION_HEADER 38 #include "multi_domain.h" 40 #include "element_with_external_element.h" 41 #include "linear_solver.h" 44 #ifdef OOMPH_HAS_TRILINOS 45 #include "trilinos_solver.h" 46 #endif // #ifdef OOMPH_HAS_TRILINOS 47 #include "iterative_linear_solver.h" 50 #include "preconditioner.h" 51 #include "general_purpose_preconditioners.h" 115 Projected_coordinate(0), Projected_lagrangian(0),
116 Projection_type(Value), Do_projection(false),
117 Backup_ninteraction(0), Backup_external_geometric_data(false) {}
129 virtual Vector<std::pair<Data*,unsigned> >
151 virtual int local_equation(
const unsigned &fld,
const unsigned &ivalue)=0;
161 (
const unsigned &fld,
const Vector<double> &s, Shape &psi)=0;
166 (
const unsigned &time,
const unsigned &fld,
const Vector<double> &s)=0;
178 template<
class ELEMENT>
181 public virtual ElementWithExternalElement
192 this->residual_for_projection
193 (residuals,GeneralisedElement::Dummy_matrix, 0);
198 ELEMENT::fill_in_contribution_to_residuals(residuals);
211 const std::string& current_string)
const 213 ElementWithExternalElement::describe_local_dofs(out,current_string);
214 ELEMENT::describe_local_dofs(out,current_string);
219 DenseMatrix<double> &jacobian)
224 this->residual_for_projection(residuals,jacobian,1);
228 ELEMENT::fill_in_contribution_to_jacobian(residuals,jacobian);
244 DenseMatrix<double> &jacobian,
245 const unsigned& flag)
248 SolidFiniteElement* solid_el_pt =
249 dynamic_cast<SolidFiniteElement*
>(
this);
251 unsigned n_dim=dim();
254 Vector<double> s(n_dim);
260 const unsigned n_node = this->nnode();
262 const unsigned n_position_type = this->nnodal_position_type();
268 const unsigned n_intpt = integral_pt()->nweight();
271 for(
unsigned ipt=0;ipt<n_intpt;ipt++)
274 for(
unsigned i=0;i<n_dim;i++) s[i] = integral_pt()->knot(ipt,i);
277 double w = integral_pt()->weight(ipt);
280 FiniteElement* other_el_pt=0;
281 other_el_pt=this->external_element_pt(0,ipt);
282 Vector<double> other_s(n_dim);
283 other_s=this->external_element_local_coord(0,ipt);
289 int local_eqn=0, local_unknown=0;
300 "Trying to project Lagrangian coordinates in non-SolidElement\n",
301 OOMPH_CURRENT_FUNCTION,
302 OOMPH_EXCEPTION_LOCATION);
306 Shape psi(n_node,n_position_type);
310 double J = this->J_eulerian(s);
318 double interpolated_xi_proj = this->interpolated_x(s,0);
321 double interpolated_xi_bar=
322 dynamic_cast<SolidFiniteElement*
>(cast_el_pt)
326 for(
unsigned l=0;l<n_node;++l)
329 for(
unsigned k=0;k<n_position_type;++k)
332 local_eqn = solid_el_pt->position_local_eqn(l,k,0);
338 residuals[local_eqn] +=
339 (interpolated_xi_proj - interpolated_xi_bar)*psi(l,k)*W;
344 for(
unsigned l2=0;l2<n_node;l2++)
347 for(
unsigned k2=0;k2<n_position_type;k2++)
350 solid_el_pt->position_local_eqn(l2,k2,0);
351 if(local_unknown >= 0)
353 jacobian(local_eqn,local_unknown)
354 += psi(l2,k2)*psi(l,k)*W;
370 Shape psi(n_node,n_position_type);
374 double J = this->J_eulerian(s);
381 double interpolated_x_proj = 0.0;
390 interpolated_x_proj = this->
get_field(0,fld,s);
394 double interpolated_x_bar=
399 for(
unsigned l=0;l<n_node;++l)
402 for(
unsigned k=0;k<n_position_type;++k)
418 if(n_position_type!=1)
421 "Trying to project generalised positions not in SolidElement\n",
422 OOMPH_CURRENT_FUNCTION,
423 OOMPH_EXCEPTION_LOCATION);
432 residuals[local_eqn] +=
433 (interpolated_x_proj - interpolated_x_bar)*psi(l,k)*W;
438 for(
unsigned l2=0;l2<n_node;l2++)
441 for(
unsigned k2=0;k2<n_position_type;k2++)
446 local_unknown = solid_el_pt
454 if(local_unknown >= 0)
456 jacobian(local_eqn,local_unknown)
457 += psi(l2,k2)*psi(l,k)*W;
483 double interpolated_value_proj = this->
get_field(now,fld,s);
486 double interpolated_value_bar =
490 for(
unsigned l=0;l<n_value;l++)
496 residuals[local_eqn] +=
497 (interpolated_value_proj - interpolated_value_bar)*psi[l]*W;
502 for(
unsigned l2=0;l2<n_value;l2++)
505 if(local_unknown >= 0)
507 jacobian(local_eqn,local_unknown)
519 "Wrong type specificied in Projection_type. This should never happen\n",
520 OOMPH_CURRENT_FUNCTION,
521 OOMPH_EXCEPTION_LOCATION);
532 const unsigned &i)
const 536 return nodal_position_gen(n,k,i);
540 return ELEMENT::zeta_nodal(n,k,i);
554 if (add_external_geometric_data())
564 if (add_external_interaction_data())
575 ignore_external_geometric_data();
576 ignore_external_interaction_data();
591 include_external_geometric_data();
595 ignore_external_geometric_data();
601 include_external_interaction_data();
605 ignore_external_interaction_data();
651 template<
class ELEMENT>
653 :
public virtual FaceGeometry<ELEMENT>
679 template<
class PROJECTABLE_ELEMENT>
686 template<
class FRIEND_PROJECTABLE_ELEMENT>
friend class RefineableTriangleMesh;
687 template<
class FRIEND_PROJECTABLE_ELEMENT>
friend class RefineableTetgenMesh;
688 template<
class FRIEND_PROJECTABLE_ELEMENT>
friend class BackupMeshForProjection;
689 template<
class FRIEND_PROJECTABLE_ELEMENT>
friend class RefineableGmshTetMesh;
703 Output_during_projection_suppressed=
true;
709 Output_during_projection_suppressed=
false;
715 {
return Use_iterative_solver_for_projection;}
719 {Use_iterative_solver_for_projection=
true;}
723 {Use_iterative_solver_for_projection=
false;}
726 void project(Mesh* base_mesh_pt,
const bool& dont_project_positions=
false)
729 if (Use_iterative_solver_for_projection)
733 #ifdef OOMPH_HAS_TRILINOS 735 if (MPI_Helpers::mpi_has_been_initialised())
738 Iterative_solver_projection_pt =
new TrilinosAztecOOSolver;
740 dynamic_cast<TrilinosAztecOOSolver*
>(Iterative_solver_projection_pt)->
741 solver_type() = TrilinosAztecOOSolver::CG;
746 Iterative_solver_projection_pt =
new CG<CRDoubleMatrix>;
750 Preconditioner_projection_pt =
new MatrixBasedDiagPreconditioner();
752 Iterative_solver_projection_pt->preconditioner_pt() =
753 Preconditioner_projection_pt;
756 Problem::linear_solver_pt() = Iterative_solver_projection_pt;
760 if (!(MPI_Helpers::mpi_has_been_initialised()))
771 Iterative_solver_projection_pt =
new CG<CRDoubleMatrix>;
774 Preconditioner_projection_pt =
new MatrixBasedDiagPreconditioner();
776 Iterative_solver_projection_pt->preconditioner_pt() =
777 Preconditioner_projection_pt;
780 Problem::linear_solver_pt() = Iterative_solver_projection_pt;
792 bool shut_up_in_newton_solve_backup=Shut_up_in_newton_solve;
795 bool backed_up_doc_time_enabled=linear_solver_pt()->is_doc_time_enabled();
796 if (Output_during_projection_suppressed)
798 linear_solver_pt()->disable_doc_time();
802 unsigned n_element = Problem::mesh_pt()->nelement();
803 unsigned n_element1=base_mesh_pt->nelement();
804 unsigned n_node = Problem::mesh_pt()->nnode();
805 unsigned n_node1=base_mesh_pt->nnode();
806 if (!Output_during_projection_suppressed)
808 oomph_info <<
"\n=============================\n";
809 oomph_info <<
"Base mesh has " << n_element1 <<
" elements\n";
810 oomph_info <<
"Target mesh has " << n_element <<
" elements\n";
811 oomph_info <<
"Base mesh has " << n_node1 <<
" nodes\n";
812 oomph_info <<
"Target mesh has " << n_node <<
" nodes\n";
813 oomph_info <<
"=============================\n\n";
819 disable_info_in_newton_solve();
826 <<
"Very odd -- no elements in target mesh; " 827 <<
" not doing anything in ProjectionProblem::project()\n";
832 unsigned nnod=Problem::mesh_pt()->nnode();
835 std::ostringstream error_stream;
837 <<
"Mesh has no nodes! Please populate the Node_pt vector\n" 838 <<
"otherwise projection won't work!\n";
839 throw OomphLibError(error_stream.str(),
840 OOMPH_CURRENT_FUNCTION,
841 OOMPH_EXCEPTION_LOCATION);
846 unsigned n_fields =
dynamic_cast<PROJECTABLE_ELEMENT*
> 850 unsigned n_dim = Problem::mesh_pt()->node_pt(0)->ndim();
853 unsigned n_history_values=0;
856 for(
unsigned e=0;e<n_element;e++)
858 PROJECTABLE_ELEMENT * el_pt =
859 dynamic_cast<PROJECTABLE_ELEMENT*
> 860 (Problem::mesh_pt()->element_pt(e));
863 el_pt->enable_projection();
869 for(
unsigned e=0;e<n_element1;e++)
871 PROJECTABLE_ELEMENT * el_pt =
872 dynamic_cast<PROJECTABLE_ELEMENT*
> 873 (base_mesh_pt->element_pt(e));
876 el_pt->enable_projection();
886 double t_start=TimingHelpers::timer();
887 Multi_domain_functions::setup_multi_domain_interaction
888 <PROJECTABLE_ELEMENT>(
this,Problem::mesh_pt(),base_mesh_pt);
889 if (!Output_during_projection_suppressed)
891 oomph_info <<
"CPU for setup of multi-domain interaction for projection: " 892 << TimingHelpers::timer()-t_start << std::endl;
894 t_start=TimingHelpers::timer();
901 if (!dont_project_positions)
910 if(dynamic_cast<SolidFiniteElement*>(
911 Problem::mesh_pt()->element_pt(0)))
914 this->store_positions();
918 this->set_coordinate_for_projection(0);
919 this->unpin_dofs_of_coordinate(0);
922 const unsigned n_lagrangian =
923 dynamic_cast<SolidNode*
>(Problem::mesh_pt()->node_pt(0))->nlagrangian();
926 for(
unsigned i=0;i<n_lagrangian;++i)
929 if (!Output_during_projection_suppressed)
931 oomph_info <<
"\n\n=============================================\n";
932 oomph_info <<
"Projecting values for Lagrangian coordinate " << i
934 oomph_info <<
"=============================================\n\n";
938 this->set_lagrangian_coordinate_for_projection(i);
941 unsigned ndof_tmp=assign_eqn_numbers();
942 if (!Output_during_projection_suppressed)
945 "Number of equations for projection of Lagrangian coordinate " 946 <<
" : "<< ndof_tmp <<std::endl << std::endl;
950 if (Problem_is_nonlinear)
952 std::ostringstream error_stream;
954 <<
"Solid mechanics problems will break if Problem_is_nonlinear is\n" 955 <<
"set to true in the projection problem because we're using the\n " 956 <<
"actual nodal positions to store the values of the Lagrangian\n" 957 <<
"coords. There shouldn't be any need for \n" 958 <<
"Problem_is_nonlinear=true anyway, apart from debugging in \n" 959 <<
"which case you now know why this case will break!\n";
960 OomphLibWarning(error_stream.str(),
961 OOMPH_CURRENT_FUNCTION,
962 OOMPH_EXCEPTION_LOCATION);
967 Problem::newton_solve();
970 unsigned n_node = Problem::mesh_pt()->nnode();
971 for(
unsigned n=0;n<n_node;++n)
974 SolidNode* solid_node_pt =
975 dynamic_cast<SolidNode*
>(Problem::mesh_pt()->node_pt(n));
977 const unsigned n_position_type = solid_node_pt->nposition_type();
979 const unsigned n_lagrangian_type = solid_node_pt->nlagrangian_type();
982 if(n_position_type != n_lagrangian_type)
984 std::ostringstream error_stream;
986 <<
"The number of generalised position dofs " 988 <<
"\n not the same as the number of generalised lagrangian dofs " 989 << n_lagrangian_type <<
".\n" 990 <<
"Projection cannot be done at the moment, sorry.\n";
992 throw OomphLibError(error_stream.str(),
993 OOMPH_CURRENT_FUNCTION,
994 OOMPH_EXCEPTION_LOCATION);
999 for(
unsigned k=0;k<n_position_type;++k)
1001 solid_node_pt->xi_gen(k,i) = solid_node_pt->x_gen(k,0);
1004 solid_node_pt->x_gen(k,0) = Solid_backup[n](k,0);
1010 this->pin_dofs_of_coordinate(0);
1015 n_history_values =
dynamic_cast<PROJECTABLE_ELEMENT*
> 1016 (Problem::mesh_pt()->element_pt(0))->
1020 if(n_history_values > 1)
1022 for (
unsigned i=0;i<n_dim;i++)
1024 if (!Output_during_projection_suppressed)
1026 oomph_info <<
"\n\n=============================================\n";
1027 oomph_info <<
"Projecting history values for coordinate " << i
1029 oomph_info <<
"=============================================\n\n";
1033 this->set_coordinate_for_projection(i);
1034 this->unpin_dofs_of_coordinate(i);
1038 for (
unsigned h_tim=n_history_values;h_tim>1;h_tim--)
1040 unsigned time_level=h_tim-1;
1043 this->set_time_level_for_projection(time_level);
1046 unsigned ndof_tmp=assign_eqn_numbers();
1047 if (!Output_during_projection_suppressed)
1049 oomph_info <<
"Number of equations for projection of coordinate " 1050 << i <<
" at time level " << time_level
1051 <<
" : "<< ndof_tmp <<std::endl << std::endl;
1055 Problem::newton_solve();
1058 unsigned n_node = Problem::mesh_pt()->nnode();
1059 for(
unsigned n=0;n<n_node;++n)
1062 Node* nod_pt = Problem::mesh_pt()->node_pt(n);
1064 const unsigned n_position_type = nod_pt->nposition_type();
1066 for(
unsigned k=0;k<n_position_type;++k)
1068 nod_pt->x_gen(time_level,k,i) = nod_pt->x_gen(0,k,i);
1071 nod_pt->x_gen(0,k,i) = Solid_backup[n](k,i);
1076 this->pin_dofs_of_coordinate(i);
1087 this->set_current_field_for_projection(0);
1088 this->unpin_dofs_of_field(0);
1099 unsigned n_element = Problem::mesh_pt()->nelement();
1100 for(
unsigned e=0;e<n_element;e++)
1102 FiniteElement* el_pt=Problem::mesh_pt()->finite_element_pt(e);
1103 unsigned nnod=el_pt->nnode();
1104 for (
unsigned j=0;j<nnod;j++)
1106 Node* nod_pt=el_pt->node_pt(j);
1107 if (nod_pt->nvalue()==0)
1109 std::ostringstream error_stream;
1110 error_stream <<
"Node at ";
1111 unsigned n=nod_pt->ndim();
1112 for (
unsigned i=0;i<n;i++)
1114 error_stream << nod_pt->x(i) <<
" ";
1117 <<
"\nhas no values. Projection (of coordinates) doesn't work\n" 1118 <<
"for such cases (at the moment), sorry! Send us the code\n" 1119 <<
"where the problem arises and we'll try to implement this\n" 1120 <<
"(up to now unnecessary) capability...\n";
1121 throw OomphLibError(error_stream.str(),
1122 OOMPH_CURRENT_FUNCTION,
1123 OOMPH_EXCEPTION_LOCATION);
1131 n_history_values =
dynamic_cast<PROJECTABLE_ELEMENT*
> 1132 (Problem::mesh_pt()->element_pt(0))->
1136 if(n_history_values > 1)
1138 for (
unsigned i=0;i<n_dim;i++)
1140 if (!Output_during_projection_suppressed)
1142 oomph_info <<
"\n\n=============================================\n";
1143 oomph_info <<
"Projecting history values for coordinate " << i
1145 oomph_info <<
"=============================================\n\n";
1149 this->set_coordinate_for_projection(i);
1153 for (
unsigned h_tim=n_history_values;h_tim>1;h_tim--)
1155 unsigned time_level=h_tim-1;
1158 this->set_time_level_for_projection(time_level);
1161 unsigned ndof_tmp=assign_eqn_numbers();
1162 if (!Output_during_projection_suppressed)
1164 oomph_info <<
"Number of equations for projection of coordinate " 1165 << i <<
" at time level " << time_level
1166 <<
" : "<< ndof_tmp <<std::endl << std::endl;
1170 Problem::newton_solve();
1173 unsigned n_element = Problem::mesh_pt()->nelement();
1174 for(
unsigned e=0;e<n_element;e++)
1176 PROJECTABLE_ELEMENT * new_el_pt =
1177 dynamic_cast<PROJECTABLE_ELEMENT*
> 1178 (Problem::mesh_pt()->element_pt(e));
1180 Vector<std::pair<Data*,unsigned> >
1181 data=new_el_pt->data_values_of_field(0);
1183 unsigned d_size=data.size();
1184 for(
unsigned d=0;d<d_size;d++)
1187 double coord=data[d].first->value(0,0);
1188 dynamic_cast<Node*
>(data[d].first)->x(time_level,i) = coord;
1196 this->pin_dofs_of_field(0);
1204 for(
unsigned e=0;e<n_element;e++)
1206 PROJECTABLE_ELEMENT * el_pt =
1207 dynamic_cast<PROJECTABLE_ELEMENT*
> 1208 (Problem::mesh_pt()->element_pt(e));
1210 el_pt->set_project_values();
1214 for (
unsigned fld=0; fld<n_fields ;fld++)
1221 this->set_current_field_for_projection(fld);
1222 this->unpin_dofs_of_field(fld);
1225 n_history_values =
dynamic_cast<PROJECTABLE_ELEMENT*
> 1230 for (
unsigned h_tim=n_history_values;h_tim>0;h_tim--)
1232 unsigned time_level=h_tim-1;
1233 if (!Output_during_projection_suppressed)
1235 oomph_info <<
"\n=========================================\n";
1236 oomph_info <<
"Projecting field " << fld <<
" at time level " 1237 << time_level<<std::endl;
1238 oomph_info <<
"========================================\n";
1242 this->set_time_level_for_projection(time_level);
1245 unsigned ndof_tmp=assign_eqn_numbers();
1246 if (!Output_during_projection_suppressed)
1248 oomph_info <<
"Number of equations for projection of field " 1249 << fld <<
" at time level " << time_level
1250 <<
" : "<< ndof_tmp <<std::endl << std::endl;
1254 Problem::newton_solve();
1261 for(
unsigned e=0;e<n_element;e++)
1263 PROJECTABLE_ELEMENT * new_el_pt =
1264 dynamic_cast<PROJECTABLE_ELEMENT*
> 1265 (Problem::mesh_pt()->element_pt(e));
1267 Vector<std::pair<Data*,unsigned> >
1268 data=new_el_pt->data_values_of_field(fld);
1270 unsigned d_size=data.size();
1271 for(
unsigned d=0;d<d_size;d++)
1274 double c_value=data[d].first->value(0,data[d].second);
1275 data[d].first->set_value(time_level,data[d].second,c_value);
1285 for(
unsigned e=0;e<n_element;e++)
1287 PROJECTABLE_ELEMENT * new_el_pt =
1288 dynamic_cast<PROJECTABLE_ELEMENT*
> 1289 (Problem::mesh_pt()->element_pt(e));
1292 new_el_pt->flush_all_external_element_storage();
1294 new_el_pt->disable_projection();
1297 for(
unsigned e=0;e<n_element1;e++)
1299 PROJECTABLE_ELEMENT * el_pt =
1300 dynamic_cast<PROJECTABLE_ELEMENT*
> 1301 (base_mesh_pt->element_pt(e));
1304 el_pt->flush_all_external_element_storage();
1307 el_pt->disable_projection();
1314 Solid_backup.clear();
1317 if (!Output_during_projection_suppressed)
1335 Shut_up_in_newton_solve=shut_up_in_newton_solve_backup;
1338 if (backed_up_doc_time_enabled)
1340 linear_solver_pt()->enable_doc_time();
1354 Problem_is_nonlinear=
false;
1358 Output_during_projection_suppressed=
true;
1361 Use_iterative_solver_for_projection=
true;
1364 Iterative_solver_projection_pt = 0;
1365 Preconditioner_projection_pt = 0;
1371 if (Iterative_solver_projection_pt!=0)
1374 delete Iterative_solver_projection_pt;
1375 Iterative_solver_projection_pt = 0;
1378 if (Preconditioner_projection_pt!=0)
1380 delete Preconditioner_projection_pt;
1381 Preconditioner_projection_pt = 0;
1393 if (Problem::mesh_pt()->nelement()==0)
return;
1398 SolidFiniteElement* solid_el_pt =
dynamic_cast<SolidFiniteElement*
> 1399 (Problem::mesh_pt()->element_pt(0));
1402 const unsigned n_node = this->mesh_pt()->nnode();
1403 Solid_backup.resize(n_node);
1405 const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1406 const unsigned n_position_type =
1407 this->mesh_pt()->node_pt(0)->nposition_type();
1409 for (
unsigned n=0;n<n_node;n++)
1412 SolidNode*
const solid_nod_pt =
1413 dynamic_cast<SolidNode*
>(this->mesh_pt()->node_pt(n));
1415 Solid_backup[n].resize(n_position_type,n_dim);
1417 for (
unsigned i=0;i<n_dim;i++)
1419 for(
unsigned k=0;k<n_position_type;k++)
1421 Solid_backup[n](k,i)= solid_nod_pt->x_gen(k,i);
1434 if (Problem::mesh_pt()->nelement()==0)
return;
1439 SolidFiniteElement* solid_el_pt =
dynamic_cast<SolidFiniteElement*
> 1440 (Problem::mesh_pt()->element_pt(0));
1443 const unsigned n_node = this->mesh_pt()->nnode();
1445 const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1446 const unsigned n_position_type =
1447 this->mesh_pt()->node_pt(0)->nposition_type();
1449 for (
unsigned n=0;n<n_node;n++)
1452 SolidNode*
const solid_nod_pt =
1453 dynamic_cast<SolidNode*
>(this->mesh_pt()->node_pt(n));
1455 for (
unsigned i=0;i<n_dim;i++)
1457 for(
unsigned k=0;k<n_position_type;k++)
1459 solid_nod_pt->x_gen(k,i) = Solid_backup[n](k,i);
1471 if (Problem::mesh_pt()->nelement()==0)
return;
1474 const unsigned n_element = Problem::mesh_pt()->nelement();
1475 for(
unsigned e=0;e<n_element;++e)
1477 FiniteElement* el_pt=Problem::mesh_pt()->finite_element_pt(e);
1478 unsigned nint=el_pt->ninternal_data();
1479 for (
unsigned j=0;j<nint;j++)
1481 Data* int_data_pt = el_pt->internal_data_pt(j);
1482 unsigned nval=int_data_pt->nvalue();
1483 for (
unsigned i=0;i<nval;i++)
1485 int_data_pt->pin(i);
1489 unsigned nnod=el_pt->nnode();
1490 for (
unsigned j=0;j<nnod;j++)
1492 Node* nod_pt = el_pt->node_pt(j);
1493 unsigned nval=nod_pt->nvalue();
1494 for (
unsigned i=0;i<nval;i++)
1502 SolidFiniteElement* solid_el_pt =
dynamic_cast<SolidFiniteElement*
> 1503 (Problem::mesh_pt()->element_pt(0));
1507 const unsigned n_node=this->mesh_pt()->nnode();
1509 if(n_node==0) {
return;}
1512 const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1513 const unsigned n_position_type =
1514 this->mesh_pt()->node_pt(0)->nposition_type();
1517 for (
unsigned n=0;n<n_node;n++)
1519 SolidNode* solid_nod_pt=
dynamic_cast<SolidNode*
>(
1520 this->mesh_pt()->node_pt(n));
1521 for (
unsigned i=0;i<n_dim;i++)
1523 for(
unsigned k=0;k<n_position_type;k++)
1525 solid_nod_pt->pin_position(k,i);
1538 if (Problem::mesh_pt()->nelement()==0)
return;
1541 const unsigned n_element = Problem::mesh_pt()->nelement();
1542 for(
unsigned e=0;e<n_element;++e)
1545 PROJECTABLE_ELEMENT * new_el_pt =
1546 dynamic_cast<PROJECTABLE_ELEMENT*
> 1547 (Problem::mesh_pt()->element_pt(e));
1549 unsigned n_fields = new_el_pt->nfields_for_projection();
1552 for(
unsigned f=0;f<n_fields;f++)
1555 Vector<std::pair<Data*,unsigned> >
1556 data=new_el_pt->data_values_of_field(f);
1558 unsigned d_size=data.size();
1559 for(
unsigned d=0;d<d_size;d++)
1561 data[d].first->unpin(data[d].second);
1567 SolidFiniteElement* solid_el_pt =
dynamic_cast<SolidFiniteElement*
> 1568 (Problem::mesh_pt()->element_pt(0));
1572 const unsigned n_node=this->mesh_pt()->nnode();
1574 if(n_node==0) {
return;}
1577 const unsigned n_dim = this->mesh_pt()->node_pt(0)->ndim();
1578 const unsigned n_position_type =
1579 this->mesh_pt()->node_pt(0)->nposition_type();
1582 for (
unsigned n=0;n<n_node;n++)
1584 SolidNode* solid_nod_pt=
dynamic_cast<SolidNode*
>(
1585 this->mesh_pt()->node_pt(n));
1586 for (
unsigned i=0;i<n_dim;i++)
1588 for(
unsigned k=0;k<n_position_type;k++)
1590 solid_nod_pt->unpin_position(k,i);
1601 const unsigned n_node = Problem::mesh_pt()->nnode();
1603 if(n_node==0) {
return;}
1606 const unsigned n_position_type =
1607 Problem::mesh_pt()->node_pt(0)->nposition_type();
1609 for(
unsigned n=0;n<n_node;++n)
1612 SolidNode* solid_nod_pt =
1613 static_cast<SolidNode*
>(Problem::mesh_pt()->node_pt(n));
1615 for(
unsigned k=0;k<n_position_type;++k)
1617 solid_nod_pt->unpin_position(k,coord);
1626 const unsigned n_node = Problem::mesh_pt()->nnode();
1628 if(n_node==0) {
return;}
1631 const unsigned n_position_type =
1632 Problem::mesh_pt()->node_pt(0)->nposition_type();
1634 for(
unsigned n=0;n<n_node;++n)
1637 SolidNode* solid_nod_pt =
1638 static_cast<SolidNode*
>(Problem::mesh_pt()->node_pt(n));
1640 for(
unsigned k=0;k<n_position_type;++k)
1642 solid_nod_pt->pin_position(k,coord);
1651 unsigned n_element = Problem::mesh_pt()->nelement();
1652 for(
unsigned e=0;e<n_element;e++)
1654 PROJECTABLE_ELEMENT * new_el_pt =
1655 dynamic_cast<PROJECTABLE_ELEMENT*
> 1656 (Problem::mesh_pt()->element_pt(e));
1658 Vector<std::pair<Data*,unsigned> >
1659 data=new_el_pt->data_values_of_field(fld);
1661 unsigned d_size=data.size();
1662 for(
unsigned d=0;d<d_size;d++)
1664 data[d].first->unpin(data[d].second);
1672 unsigned n_element = Problem::mesh_pt()->nelement();
1673 for(
unsigned e=0;e<n_element;e++)
1675 PROJECTABLE_ELEMENT * new_el_pt =
1676 dynamic_cast<PROJECTABLE_ELEMENT*
> 1677 (Problem::mesh_pt()->element_pt(e));
1679 Vector<std::pair<Data*,unsigned> >
1680 data=new_el_pt->data_values_of_field(fld);
1682 unsigned d_size=data.size();
1683 for(
unsigned d=0;d<d_size;d++)
1685 data[d].first->pin(data[d].second);
1693 unsigned n_element = Problem::mesh_pt()->nelement();
1694 for(
unsigned e=0;e<n_element;e++)
1696 PROJECTABLE_ELEMENT * el_pt =
1697 dynamic_cast<PROJECTABLE_ELEMENT*
> 1698 (Problem::mesh_pt()->element_pt(e));
1701 el_pt->time_level_for_projection()=time_level;
1708 unsigned n_element = Problem::mesh_pt()->nelement();
1709 for(
unsigned e=0;e<n_element;e++)
1711 PROJECTABLE_ELEMENT * new_el_pt =
1712 dynamic_cast<PROJECTABLE_ELEMENT*
> 1713 (Problem::mesh_pt()->element_pt(e));
1716 new_el_pt->set_project_coordinates();
1719 new_el_pt->projected_coordinate()=i;
1726 unsigned n_element = Problem::mesh_pt()->nelement();
1727 for(
unsigned e=0;e<n_element;e++)
1729 PROJECTABLE_ELEMENT * new_el_pt =
1730 dynamic_cast<PROJECTABLE_ELEMENT*
> 1731 (Problem::mesh_pt()->element_pt(e));
1734 new_el_pt->set_project_lagrangian();
1737 new_el_pt->projected_lagrangian_coordinate()=i;
1744 unsigned n_element = Problem::mesh_pt()->nelement();
1745 for(
unsigned e=0;e<n_element;e++)
1747 PROJECTABLE_ELEMENT * new_el_pt =
1748 dynamic_cast<PROJECTABLE_ELEMENT*
> 1749 (Problem::mesh_pt()->element_pt(e));
1752 new_el_pt->projected_field()=fld;
unsigned Backup_ninteraction
Store number of "external" interactions that were assigned to the element before doing the projection...
Projection_Type Projection_type
Variable to indicate if we're projecting the history values of the nodal coordinates (Coordinate) the...
void set_coordinate_for_projection(const unsigned &i)
Set the coordinate for projection.
void disable_projection()
Helper function to restore the element to the state it was in before we entered the projection mode a...
void unpin_dofs_of_coordinate(const unsigned &coord)
Helper function to unpin the values of coordinate coord.
unsigned & time_level_for_projection()
Which history value are we projecting?
IterativeLinearSolver * Iterative_solver_projection_pt
void pin_all()
Pin all the field values and position unknowns (bit inefficient)
ProjectableElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
virtual double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)=0
Return Jacobian of mapping and the shape functions associated with field fld. The number of shape fun...
void set_project_lagrangian()
Project (current and only values of) lagrangian coordinates.
bool Do_projection
Bool to know if we do projection or not. If false (the default) we solve the element's "real" equatio...
void pin_dofs_of_field(const unsigned &fld)
Helper function to unpin dofs of fld-th field.
Template-free Base class for projectable elements.
void enable_suppress_output_during_projection()
Suppress all output during projection phases.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Overloaded version of fill_in_contribution_to_residuals.
unsigned & projected_field()
Field that is currently being projected.
void set_project_coordinates()
Project (history values of) coordintes.
void project(Mesh *base_mesh_pt, const bool &dont_project_positions=false)
Project from base into the problem's own mesh.
unsigned Time_level_for_projection
Time level we are projecting (0=current values; >0: history values)
virtual ~ProjectableElementBase()
Virtual destructor.
void set_lagrangian_coordinate_for_projection(const unsigned &i)
Set the Lagrangian coordinate for projection.
void store_positions()
Helper function to store positions (the only things that have been set before doing projection...
unsigned Projected_coordinate
When projecting the history values of the nodal coordinates, this is the coordinate we're projecting...
bool Output_during_projection_suppressed
Flag to suppress output during projection.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Overloaded version of fill_in_contribution_to_jacobian.
void unpin_dofs_of_field(const unsigned &fld)
Helper function to unpin dofs of fld-th field.
bool use_iterative_solver_for_projection()
virtual void residual_for_projection(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Residual for the projection step. Flag indicates if we want the Jacobian (1) or not (0)...
virtual unsigned nhistory_values_for_projection(const unsigned &fld)=0
Number of history values to be stored for fld-th field (includes current value!)
virtual unsigned nvalue_of_field(const unsigned &fld)=0
Return number of values (pinned or not) associated with field fld within the element. This must correspond to the number of shape functions returned in jacobian_and_shape_of_field(...).
virtual unsigned nhistory_values_for_coordinate_projection()=0
Number of history values to be stored when projecting the history values of the nodal coordinates (in...
void enable_projection()
Backup the element's state and switch it to projection mode.
bool Backup_external_interaction_data
Remember if the element includes external data when used in non-projection mode (this is temporarily ...
void enable_use_iterative_solver_for_projection()
Enables the use of an iterative solver for projection.
unsigned & projected_coordinate()
When projecting the history values of the nodal coordinates, this is the coordinate we're projecting...
void pin_dofs_of_coordinate(const unsigned &coord)
Helper function to unpin the values of coordinate coord.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Use Eulerian coordinates for matching in locate_zeta when doing projection.
unsigned & projected_lagrangian_coordinate()
When projecting the Lagrangian coordinates this is the coordinate that is being projected.
void unpin_all()
Unpin all the field values and position unknowns (bit inefficient)
virtual unsigned nfields_for_projection()=0
Number of fields of the problem, so e.g. for 2D Navier Stokes this would be 3 (for the two velocities...
void set_time_level_for_projection(const unsigned &time_level)
Helper function to set time level for projection.
void disable_use_iterative_solver_for_projection()
Disbales the use of an iterative solver for projection.
void set_current_field_for_projection(const unsigned &fld)
Set current field for projection.
Vector< DenseMatrix< double > > Solid_backup
Backup for pin status of solid node's position Data.
virtual double get_field(const unsigned &time, const unsigned &fld, const Vector< double > &s)=0
Return the fld-th field at local coordinates s at time-level time (time=0: current value; time>0: his...
void disable_suppress_output_during_projection()
Undo suppression of all output during projection phases.
void restore_positions()
Helper function to restore positions (the only things that have been set before doing projection...
virtual Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)=0
Pure virtual function in which the element writer must specify the values associated with field fld...
unsigned Projected_field
Field that is currently being projected.
virtual int local_equation(const unsigned &fld, const unsigned &ivalue)=0
Return local equation numbers associated with value ivalue of field fld within the element...
bool Use_iterative_solver_for_projection
void describe_local_dofs(std::ostream &out, const std::string ¤t_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
Projection_Type
Enumerated collection to specify which projection problem is to be solved.
bool Backup_external_geometric_data
Remember if the element includes external geometric data when used in non-projection mode (this is te...
Preconditioner * Preconditioner_projection_pt
void set_project_values()
Project (history values of) values.
unsigned Projected_lagrangian
When projecting the Lagrangain coordinates indicate which coordiante is to be projected.