78 for (
unsigned i = 0;
i < n_internal_data;
i++)
88 <<
"---------------------------------------------------------------"
91 <<
"Info: Data is already included in this element's internal "
94 <<
"It's stored as entry " <<
i <<
" and I'm not adding it again."
97 <<
"Note: You can suppress this message by recompiling without"
98 <<
"\n PARANOID or setting the boolean \n"
100 "GeneralisedElement::Suppress_warning_about_repeated_internal_"
102 <<
"\n\n to true." << std::endl
103 <<
"---------------------------------------------------------------"
115 Data** new_data_pt =
new Data*[n_internal_data + n_external_data + 1];
118 for (
unsigned i = 0;
i < n_internal_data;
i++)
124 new_data_pt[n_internal_data] = data_pt;
127 for (
unsigned i = 0;
i < n_external_data;
i++)
129 new_data_pt[n_internal_data + 1 +
i] =
Data_pt[n_internal_data +
i];
139 Data_fd.resize(n_internal_data + n_external_data + 1);
141 for (
unsigned i = n_external_data;
i > 0;
i--)
152 return n_internal_data;
162 std::deque<unsigned long>
const& global_eqn_numbers,
163 std::deque<double*>
const& global_dof_pt)
166 const unsigned n_dof =
Ndof;
168 const unsigned n_additional_dof = global_eqn_numbers.size();
170 if (n_additional_dof == 0)
176 const unsigned new_n_dof = n_dof + n_additional_dof;
178 unsigned long* new_eqn_number =
new unsigned long[new_n_dof];
181 for (
unsigned i = 0;
i < n_dof;
i++)
187 unsigned index = n_dof;
189 for (std::deque<unsigned long>::const_iterator it =
190 global_eqn_numbers.begin();
191 it != global_eqn_numbers.end();
195 new_eqn_number[index] = *it;
202 const unsigned n_additional_dof_pt = global_dof_pt.size();
203 if (n_additional_dof_pt > 0)
207 if (n_additional_dof_pt != n_additional_dof)
209 std::ostringstream error_stream;
211 <<
"global_dof_pt is non-empty, yet it does not have the same size\n"
212 <<
"as global_eqn_numbers.\n"
213 <<
"There are " << n_additional_dof <<
" equation numbers,\n"
214 <<
"but " << n_additional_dof_pt << std::endl;
217 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
222 double** new_dof_pt =
new double*[new_n_dof];
224 for (
unsigned i = 0;
i < n_dof;
i++)
230 unsigned index = n_dof;
232 for (std::deque<double*>::const_iterator it = global_dof_pt.begin();
233 it != global_dof_pt.end();
237 new_dof_pt[index] = *it;
281 GeneralisedElement::~GeneralisedElement()
322 for (
unsigned i = 0;
i < n_external_data;
i++)
332 <<
"---------------------------------------------------------------"
335 <<
"Info: Data is already included in this element's external "
338 <<
"It's stored as entry " <<
i <<
" and I'm not adding it again"
341 <<
"Note: You can suppress this message by recompiling without"
342 <<
"\n PARANOID or setting the boolean \n"
344 "GeneralisedElement::Suppress_warning_about_repeated_external_"
346 <<
"\n\n to true." << std::endl
347 <<
"---------------------------------------------------------------"
359 Data** new_data_pt =
new Data*[n_internal_data + n_external_data + 1];
362 for (
unsigned i = 0;
i < (n_internal_data + n_external_data);
i++)
368 new_data_pt[n_internal_data + n_external_data] = data_pt;
377 Data_fd.resize(n_internal_data + n_external_data + 1);
379 Data_fd[n_internal_data + n_external_data] = fd;
385 return n_external_data;
397 if (n_external_data > 0)
400 Data** new_data_pt = 0;
406 if (n_internal_data > 0)
409 new_data_pt =
new Data*[n_internal_data];
411 for (
unsigned i = 0;
i < n_internal_data;
i++)
425 Data_fd.resize(n_internal_data);
441 unsigned index = n_external_data;
443 for (
unsigned i = 0;
i < n_external_data;
i++)
455 if (index < n_external_data)
458 Data** new_data_pt = 0;
464 const unsigned n_total_data = n_internal_data + n_external_data - 1;
467 if (n_total_data > 0)
469 new_data_pt =
new Data*[n_total_data];
473 for (
unsigned i = 0;
i < n_internal_data;
i++)
479 unsigned counter = 0;
480 for (
unsigned i = 0;
i < n_external_data;
i++)
487 new_data_pt[n_internal_data + counter] =
Data_pt[
i];
508 std::ostringstream warning_stream;
510 <<
"Data removed from element's external data " << std::endl
511 <<
"You may have to update the indices for remaining data "
513 <<
"This can be achieved by using add_external_data() "
516 "GeneralisedElement::flush_external_data()",
517 OOMPH_EXCEPTION_LOCATION);
557 std::ostream& out,
const std::string& current_string)
const
561 std::stringstream conversion;
562 conversion <<
" of Internal Data " <<
i << current_string;
579 std::ostream& out,
const std::string& current_string)
const
586 for (
unsigned i = 0;
i < n_internal_data;
i++)
591 std::stringstream conversion;
592 conversion <<
" of Internal Data " <<
i << current_string;
599 for (
unsigned i = 0;
i < n_external_data;
i++)
604 std::stringstream conversion;
605 conversion <<
" of External Data " <<
i << current_string;
617 std::map<unsigned, double*>& map_of_value_pt)
676 const Vector<long>& vector_of_eqn_numbers,
unsigned& index)
697 const bool& store_local_dof_pt)
706 std::ostringstream error_stream;
710 std::map<unsigned, bool> is_repeated;
711 std::set<unsigned long> set_of_global_eqn_numbers;
712 unsigned old_ndof = 0;
713 for (
unsigned n = 0; n <
Ndof; ++n)
715 set_of_global_eqn_numbers.insert(
Eqn_number[n]);
716 if (set_of_global_eqn_numbers.size() == old_ndof)
718 error_stream <<
"Repeated global eqn: " <<
Eqn_number[n] << std::endl;
721 old_ndof = set_of_global_eqn_numbers.size();
725 if (set_of_global_eqn_numbers.size() !=
Ndof)
728 error_stream <<
"Element is ";
729 if (!
is_halo()) error_stream <<
"not ";
730 error_stream <<
"a halo element\n\n";
732 error_stream <<
"\nLocal/lobal equation numbers: " << std::endl;
733 for (
unsigned n = 0; n <
Ndof; ++n)
735 error_stream <<
" " << n <<
" " <<
Eqn_number[n] << std::endl;
739 error_stream << std::endl <<
" element address is " <<
this << std::endl;
743 error_stream <<
"Number of internal data " << nint << std::endl;
744 for (
unsigned i = 0;
i < nint;
i++)
747 unsigned nval = data_pt->
nvalue();
748 for (
unsigned j = 0; j < nval; j++)
751 error_stream <<
"Internal dof: " << eqn_no << std::endl;
752 if (is_repeated[
unsigned(eqn_no)])
754 error_stream <<
"Repeated internal dof: " << eqn_no << std::endl;
761 error_stream <<
"Number of external data " << next << std::endl;
762 for (
unsigned i = 0;
i < next;
i++)
765 unsigned nval = data_pt->
nvalue();
766 for (
unsigned j = 0; j < nval; j++)
769 error_stream <<
"External dof: " << eqn_no << std::endl;
770 if (is_repeated[
unsigned(eqn_no)])
772 error_stream <<
"Repeated external dof: " << eqn_no;
773 Node* nod_pt =
dynamic_cast<Node*
>(data_pt);
776 error_stream <<
" (is a node at: ";
777 unsigned ndim = nod_pt->
ndim();
778 for (
unsigned ii = 0; ii < ndim; ii++)
780 error_stream << nod_pt->
x(ii) <<
" ";
783 error_stream <<
")\n";
798 error_stream <<
"Number of external element data " << next
801 for (
unsigned i = 0;
i < next;
i++)
803 unsigned nval = data_pt[
i]->nvalue();
804 for (
unsigned j = 0; j < nval; j++)
806 int eqn_no = data_pt[
i]->eqn_number(j);
807 error_stream <<
"External element dof: " << eqn_no << std::endl;
808 if (is_repeated[
unsigned(eqn_no)])
810 error_stream <<
"Repeated external element dof: " << eqn_no;
811 Node* nod_pt =
dynamic_cast<Node*
>(data_pt[
i]);
814 error_stream <<
" (is a node at: ";
815 unsigned ndim = nod_pt->
ndim();
816 for (
unsigned ii = 0; ii < ndim; ii++)
818 error_stream << nod_pt->
x(ii) <<
" ";
821 error_stream <<
")\n";
831 error_stream <<
"Number of external element geom data " << next
835 for (
unsigned i = 0;
i < next;
i++)
837 unsigned nval = data_pt[
i]->nvalue();
838 for (
unsigned j = 0; j < nval; j++)
840 int eqn_no = data_pt[
i]->eqn_number(j);
841 error_stream <<
"External element geometric dof: " << eqn_no
843 if (is_repeated[
unsigned(eqn_no)])
846 <<
"Repeated external element geometric dof: " << eqn_no
847 <<
" " << data_pt[
i]->value(j);
848 Node* nod_pt =
dynamic_cast<Node*
>(data_pt[
i]);
851 error_stream <<
" (is a node at: ";
852 unsigned ndim = nod_pt->
ndim();
853 for (
unsigned ii = 0; ii < ndim; ii++)
855 error_stream << nod_pt->
x(
i) <<
" ";
859 error_stream <<
"\n";
870 unsigned n_node = f_el_pt->
nnode();
871 for (
unsigned n = 0; n < n_node; n++)
874 unsigned nval = nod_pt->
nvalue();
875 for (
unsigned j = 0; j < nval; j++)
878 error_stream <<
"Node " << n <<
": Nodal dof: " << eqn_no;
881 if (is_repeated[
unsigned(eqn_no)])
883 error_stream <<
"Node " << n
884 <<
": Repeated nodal dof: " << eqn_no;
887 error_stream <<
" (resized)";
889 error_stream << std::endl;
894 if (solid_nod_pt != 0)
897 unsigned nval = data_pt->
nvalue();
898 for (
unsigned j = 0; j < nval; j++)
901 error_stream <<
"Node " << n <<
": Positional dof: " << eqn_no
903 if (is_repeated[
unsigned(eqn_no)])
905 error_stream <<
"Repeated positional dof: " << eqn_no <<
" "
906 << data_pt->
value(j) << std::endl;
913 n_node = f_el_pt->
nnode();
914 for (
unsigned n = 0; n < n_node; n++)
917 unsigned n_dim = nod_pt->
ndim();
918 error_stream <<
"Node " << n <<
" at ( ";
919 for (
unsigned i = 0;
i < n_dim;
i++)
921 error_stream << nod_pt->
x(
i) <<
" ";
923 error_stream <<
")" << std::endl;
931 error_stream << std::endl << std::endl;
933 <<
"---------------------------------------------------------------"
936 <<
"Note: You can suppress this warning by recompiling without"
937 <<
"\n PARANOID or setting the boolean \n"
939 "GeneralisedElement::Suppress_warning_about_any_repeated_data"
940 <<
"\n\n to true." << std::endl
942 <<
"Only do this if you know what you're doing; repeated equation\n"
943 <<
"numbers are usually a sign of trouble...\n"
944 <<
"---------------------------------------------------------------"
951 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
967 const bool& store_local_dof_pt)
973 const unsigned n_total_data = n_internal_data + n_external_data;
976 if (n_total_data > 0)
986 for (
unsigned i = 1;
i < n_total_data; ++
i)
1000 if (n_total_values == 0)
1015 for (
unsigned i = 0;
i < n_total_values; ++
i)
1021 for (
unsigned i = 1;
i < n_total_data; ++
i)
1032 std::deque<unsigned long> global_eqn_number_queue;
1035 for (
unsigned i = 0;
i < n_internal_data;
i++)
1040 unsigned n_value = data_pt->
nvalue();
1043 for (
unsigned j = 0; j < n_value; j++)
1051 global_eqn_number_queue.push_back(
eqn_number);
1053 if (store_local_dof_pt)
1072 for (
unsigned i = 0;
i < n_external_data;
i++)
1077 unsigned n_value = data_pt->
nvalue();
1080 for (
unsigned j = 0; j < n_value; j++)
1088 global_eqn_number_queue.push_back(
eqn_number);
1090 if (store_local_dof_pt)
1111 if (store_local_dof_pt)
1133 const bool& fd_all_data)
1139 if (n_internal_data == 0)
1149 const unsigned n_dof =
ndof();
1155 int local_unknown = 0;
1161 for (
unsigned i = 0;
i < n_internal_data;
i++)
1171 for (
unsigned j = 0; j < n_value; j++)
1176 if (local_unknown >= 0)
1182 const double old_var = *value_pt;
1185 *value_pt += fd_step;
1194 for (
unsigned m = 0; m < n_dof; m++)
1196 double sum = (newres[m] - residuals[m]) / fd_step;
1198 jacobian(m, local_unknown) = sum;
1202 *value_pt = old_var;
1230 const bool& fd_all_data)
1235 if (n_external_data == 0)
1245 const unsigned n_dof =
ndof();
1251 int local_unknown = 0;
1257 for (
unsigned i = 0;
i < n_external_data;
i++)
1267 for (
unsigned j = 0; j < n_value; j++)
1272 if (local_unknown >= 0)
1278 const double old_var = *value_pt;
1281 *value_pt += fd_step;
1290 for (
unsigned m = 0; m < n_dof; m++)
1292 double sum = (newres[m] - residuals[m]) / fd_step;
1294 jacobian(m, local_unknown) = sum;
1298 *value_pt = old_var;
1324 "Empty fill_in_contribution_to_mass_matrix() has been called.\n";
1326 "This function is called from the default implementation of\n";
1327 error_message +=
"get_mass_matrix();\n";
1328 error_message +=
"and must calculate the residuals vector and mass matrix ";
1329 error_message +=
"without initialising any of their entries.\n\n";
1332 "If you do not wish to use these defaults, you must overload\n";
1333 error_message +=
"get_mass_matrix(), which must initialise the entries\n";
1334 error_message +=
"of the residuals vector and mass matrix to zero.\n";
1338 "GeneralisedElement::fill_in_contribution_to_mass_matrix()",
1339 OOMPH_EXCEPTION_LOCATION);
1356 "Empty fill_in_contribution_to_jacobian_and_mass_matrix() has been ";
1357 error_message +=
"called.\n";
1359 "This function is called from the default implementation of\n";
1360 error_message +=
"get_jacobian_and_mass_matrix();\n";
1362 "and must calculate the residuals vector and mass and jacobian matrices ";
1363 error_message +=
"without initialising any of their entries.\n\n";
1366 "If you do not wish to use these defaults, you must overload\n";
1368 "get_jacobian_and_mass_matrix(), which must initialise the entries\n";
1370 "of the residuals vector, jacobian and mass matrix to zero.\n";
1374 "GeneralisedElement::fill_in_contribution_to_jacobian_and_mass_matrix()",
1375 OOMPH_EXCEPTION_LOCATION);
1390 "Empty fill_in_contribution_to_dresiduals_dparameter() has been ";
1391 error_message +=
"called.\n";
1393 "This function is called from the default implementation of\n";
1394 error_message +=
"get_dresiduals_dparameter();\n";
1395 error_message +=
"and must calculate the derivatives of the residuals "
1396 "vector with respect\n";
1397 error_message +=
"to the passed parameter ";
1398 error_message +=
"without initialising any values.\n\n";
1401 "If you do not wish to use these defaults, you must overload\n";
1403 "get_dresiduals_dparameter(), which must initialise the entries\n";
1404 error_message +=
"of the returned vector to zero.\n";
1407 "This function is intended for use instead of the default (global) \n";
1409 "finite-difference routine when analytic expressions are to be used\n";
1410 error_message +=
"in continuation and bifurcation tracking problems.\n\n";
1411 error_message +=
"This function is only called when the function\n";
1413 "Problem::set_analytic_dparameter() has been called in the driver code\n";
1417 "GeneralisedElement::fill_in_contribution_to_dresiduals_dparameter()",
1418 OOMPH_EXCEPTION_LOCATION);
1431 double*
const& parameter_pt,
1436 "Empty fill_in_contribution_to_djacobian_dparameter() has been ";
1437 error_message +=
"called.\n";
1439 "This function is called from the default implementation of\n";
1440 error_message +=
"get_djacobian_dparameter();\n";
1442 "and must calculate the derivatives of residuals vector and jacobian ";
1443 error_message +=
"matrix\n";
1444 error_message +=
"with respect to the passed parameter ";
1445 error_message +=
"without initialising any values.\n\n";
1448 "If you do not wish to use these defaults, you must overload\n";
1450 "get_djacobian_dparameter(), which must initialise the entries\n";
1451 error_message +=
"of the returned vector and matrix to zero.\n\n";
1454 "This function is intended for use instead of the default (global) \n";
1456 "finite-difference routine when analytic expressions are to be used\n";
1457 error_message +=
"in continuation and bifurcation tracking problems.\n\n";
1458 error_message +=
"This function is only called when the function\n";
1460 "Problem::set_analytic_dparameter() has been called in the driver code\n";
1465 "GeneralisedElement::fill_in_contribution_to_djacobian_dparameter()",
1466 OOMPH_EXCEPTION_LOCATION);
1482 double*
const& parameter_pt,
1488 "Empty fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter() "
1490 error_message +=
" been called.\n";
1492 "This function is called from the default implementation of\n";
1493 error_message +=
"get_djacobian_and_dmass_matrix_dparameter();\n";
1495 "and must calculate the residuals vector and mass and jacobian matrices ";
1496 error_message +=
"without initialising any of their entries.\n\n";
1499 "If you do not wish to use these defaults, you must overload\n";
1500 error_message +=
"get_djacobian_and_dmass_matrix_dparameter(), which must "
1502 error_message +=
"entries of the returned vector and matrices to zero.\n";
1506 "This function is intended for use instead of the default (global) \n";
1508 "finite-difference routine when analytic expressions are to be used\n";
1509 error_message +=
"in continuation and bifurcation tracking problems.\n\n";
1510 error_message +=
"This function is only called when the function\n";
1512 "Problem::set_analytic_dparameter() has been called in the driver code\n";
1516 "GeneralisedElement::fill_in_contribution_to_djacobian_"
1517 "and_dmass_matrix_dparameter()",
1518 OOMPH_EXCEPTION_LOCATION);
1537 "Empty fill_in_contribution_to_hessian_vector_products() has been ";
1538 error_message +=
"called.\n";
1540 "This function is called from the default implementation of\n";
1541 error_message +=
"get_hessian_vector_products(); ";
1542 error_message +=
" and must calculate the products \n";
1543 error_message +=
"of the hessian matrix with the (global) ";
1544 error_message +=
"vectors Y and C\n";
1545 error_message +=
"without initialising any values.\n\n";
1548 "If you do not wish to use these defaults, you must overload\n";
1550 "get_hessian_vector_products(), which must initialise the entries\n";
1551 error_message +=
"of the returned vector to zero.\n\n";
1554 "This function is intended for use instead of the default (global) \n";
1556 "finite-difference routine when analytic expressions are to be used.\n";
1557 error_message +=
"This function is only called when the function\n";
1558 error_message +=
"Problem::set_analytic_hessian_products() has been called "
1559 "in the driver code\n";
1563 "GeneralisedElement::fill_in_contribution_to_hessian_vector_product()",
1564 OOMPH_EXCEPTION_LOCATION);
1572 Vector<std::pair<unsigned, unsigned>>
const& history_index,
1576 "Empty fill_in_contribution_to_inner_products() has been called.\n";
1578 "This function is called from the default implementation of\n";
1579 error_message +=
"get_inner_products();\n ";
1580 error_message +=
"It must calculate the inner products over the element\n";
1581 error_message +=
"of given pairs of history values\n";
1582 error_message +=
"without initialision of the output vector.\n\n";
1585 "If you do not wish to use these defaults, you must overload\n";
1587 "get_inner_products(), which must initialise the entries\n";
1588 error_message +=
"of the returned vector to zero.\n\n";
1591 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1604 "Empty fill_in_contribution_to_inner_product_vectors() has been "
1607 "This function is called from the default implementation of\n";
1608 error_message +=
"get_inner_product_vectors(); ";
1609 error_message +=
" and must calculate the vectors \n";
1610 error_message +=
"corresponding to the input history values\n ";
1612 "that when multiplied by other vectors of history values\n";
1613 error_message +=
"return the appropriate dot products.\n\n";
1614 error_message +=
"The output vectors must not be initialised.\n\n";
1617 "If you do not wish to use these defaults, you must overload\n";
1619 "get_inner_products(), which must initialise the entries\n";
1620 error_message +=
"of the returned vector to zero.\n\n";
1623 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1642 oomph_info <<
"\n ERROR: Failed GeneralisedElement::self_test()!"
1644 oomph_info <<
"for internal data object number: " <<
i << std::endl;
1654 oomph_info <<
"\n ERROR: Failed GeneralisedElement::self_test()!"
1656 oomph_info <<
"for external data object number: " <<
i << std::endl;
1676 namespace Locate_zeta_helpers
1715 const unsigned&
i)
const
1738 std::ostream& out,
const std::string& current_string)
const
1756 std::ostream& out,
const std::string& current_string)
const
1759 const unsigned n_node =
nnode();
1760 for (
unsigned n = 0; n < n_node; n++)
1765 std::stringstream conversion;
1766 conversion <<
" of Node " << n << current_string;
1789 std::ostringstream error_stream;
1791 <<
"Determinant of Jacobian matrix is zero --- "
1792 <<
"singular mapping!\nThe determinant of the "
1793 <<
"jacobian is " << std::fabs(jacobian)
1794 <<
" which is smaller than the treshold "
1796 <<
"You can change this treshold, by specifying "
1797 <<
"FiniteElement::Tolerance_for_singular_jacobian \n"
1798 <<
"Here are the nodal coordinates of the inverted element\n"
1799 <<
"in the form \n\n x,y[,z], hang_status\n\n"
1800 <<
"where hang_status = 1 or 2 for non-hanging or hanging\n"
1801 <<
"nodes respectively (useful for automatic sizing of\n"
1802 <<
"tecplot markers to identify the hanging nodes). \n\n";
1804 unsigned n_nod =
nnode();
1805 unsigned hang_count = 0;
1806 for (
unsigned j = 0; j < n_nod; j++)
1808 for (
unsigned i = 0;
i < n_dim_nod;
i++)
1810 error_stream <<
node_pt(j)->
x(
i) <<
" ";
1814 error_stream <<
" 2";
1819 error_stream <<
" 1";
1821 error_stream << std::endl;
1823 error_stream << std::endl << std::endl;
1827 <<
"NOTE: Offending element is associated with a MacroElement\n"
1828 <<
" AND the element has hanging nodes! \n"
1829 <<
" If an element is thin and highly curved, the \n"
1830 <<
" constraints imposed by\n \n"
1831 <<
" (1) inter-element continuity (imposed by the hanging\n"
1832 <<
" node constraints) and \n\n"
1833 <<
" (2) the need to respect curvilinear domain boundaries\n"
1834 <<
" during mesh refinement (imposed by the element's\n"
1835 <<
" macro element mapping)\n\n"
1836 <<
" may be irreconcilable! \n \n"
1837 <<
" You may have to re-design your base mesh to avoid \n"
1838 <<
" the creation of thin, highly curved elements during\n"
1839 <<
" the refinement process.\n"
1843 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1857 std::ostringstream error_stream;
1858 error_stream <<
"Negative Jacobian in transform from "
1859 <<
"local to global coordinates" << std::endl
1860 <<
" You have an inverted coordinate system!"
1864 <<
"Here are the nodal coordinates of the inverted element\n"
1865 <<
"in the form \n\n x,y[,z], hang_status\n\n"
1866 <<
"where hang_status = 1 or 2 for non-hanging or hanging\n"
1867 <<
"nodes respectively (useful for automatic sizing of\n"
1868 <<
"tecplot markers to identify the hanging nodes). \n\n";
1870 unsigned n_nod =
nnode();
1871 unsigned hang_count = 0;
1872 for (
unsigned j = 0; j < n_nod; j++)
1874 for (
unsigned i = 0;
i < n_dim_nod;
i++)
1876 error_stream <<
node_pt(j)->
x(
i) <<
" ";
1880 error_stream <<
" 2";
1885 error_stream <<
" 1";
1887 error_stream << std::endl;
1889 error_stream << std::endl << std::endl;
1893 <<
"NOTE: The inverted element is associated with a MacroElement\n"
1894 <<
" AND the element has hanging nodes! \n"
1895 <<
" If an element is thin and highly curved, the \n"
1896 <<
" constraints imposed by\n \n"
1897 <<
" (1) inter-element continuity (imposed by the hanging\n"
1898 <<
" node constraints) and \n\n"
1899 <<
" (2) the need to respect curvilinear domain boundaries\n"
1900 <<
" during mesh refinement (imposed by the element's\n"
1901 <<
" macro element mapping)\n\n"
1902 <<
" may be irreconcilable! \n \n"
1903 <<
" You may have to re-design your base mesh to avoid \n"
1904 <<
" the creation of thin, highly curved elements during\n"
1905 <<
" the refinement process.\n"
1912 <<
"If you believe that inverted elements do not cause any\n"
1913 <<
"problems in your specific application you can \n "
1914 <<
"suppress this test by: " << std::endl
1915 <<
" i) setting the (static) flag "
1916 <<
"FiniteElement::Accept_negative_jacobian to be true" << std::endl;
1917 error_stream <<
" ii) switching OFF the PARANOID flag" << std::endl
1923 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1939 const unsigned el_dim =
dim();
1942 const unsigned n_shape =
nnode();
1951 std::ostringstream error_message;
1952 error_message <<
"Dimension mismatch" << std::endl;
1955 <<
" for the jacobian of the mapping to be well-defined"
1958 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1964 for (
unsigned i = 0;
i < el_dim;
i++)
1967 for (
unsigned j = 0; j < el_dim; j++)
1970 jacobian(
i, j) = 0.0;
1972 for (
unsigned l = 0; l < n_shape; l++)
1974 for (
unsigned k = 0; k < n_shape_type; k++)
1997 const unsigned el_dim =
dim();
2001 const unsigned n_shape =
nnode();
2004 const unsigned n_row =
N2deriv[el_dim];
2009 for (
unsigned i = 0;
i < n_row;
i++)
2012 for (
unsigned j = 0; j < el_dim; j++)
2015 jacobian2(
i, j) = 0.0;
2017 for (
unsigned l = 0; l < n_shape; l++)
2020 for (
unsigned k = 0; k < n_shape_type; k++)
2043 const unsigned n_node =
nnode();
2047 const unsigned n_dim_element =
dim();
2051 for (
unsigned i = 0;
i < n_dim_element;
i++)
2053 for (
unsigned j = 0; j < n_dim_node; j++)
2056 interpolated_G(
i, j) = 0.0;
2057 for (
unsigned l = 0; l < n_node; l++)
2059 for (
unsigned k = 0; k < n_position_type; k++)
2061 interpolated_G(
i, j) +=
2074 double FiniteElement::invert_jacobian<0>(
2079 oomph_info <<
"\nWarning: You are trying to invert the jacobian for "
2080 <<
"a 'point' element" << std::endl
2081 <<
"This makes no sense and is almost certainly an error"
2094 double FiniteElement::invert_jacobian<1>(
2099 const double det = jacobian(0, 0);
2107 inverse_jacobian(0, 0) = 1.0 / jacobian(0, 0);
2117 double FiniteElement::invert_jacobian<2>(
2123 jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0);
2131 inverse_jacobian(0, 0) = jacobian(1, 1) / det;
2132 inverse_jacobian(0, 1) = -jacobian(0, 1) / det;
2133 inverse_jacobian(1, 0) = -jacobian(1, 0) / det;
2134 inverse_jacobian(1, 1) = jacobian(0, 0) / det;
2145 double FiniteElement::invert_jacobian<3>(
2150 const double det = jacobian(0, 0) * jacobian(1, 1) * jacobian(2, 2) +
2151 jacobian(0, 1) * jacobian(1, 2) * jacobian(2, 0) +
2152 jacobian(0, 2) * jacobian(1, 0) * jacobian(2, 1) -
2153 jacobian(0, 0) * jacobian(1, 2) * jacobian(2, 1) -
2154 jacobian(0, 1) * jacobian(1, 0) * jacobian(2, 2) -
2155 jacobian(0, 2) * jacobian(1, 1) * jacobian(2, 0);
2163 inverse_jacobian(0, 0) =
2164 (jacobian(1, 1) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 1)) / det;
2165 inverse_jacobian(0, 1) =
2166 -(jacobian(0, 1) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 1)) /
2168 inverse_jacobian(0, 2) =
2169 (jacobian(0, 1) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 1)) / det;
2170 inverse_jacobian(1, 0) =
2171 -(jacobian(1, 0) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 0)) /
2173 inverse_jacobian(1, 1) =
2174 (jacobian(0, 0) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 0)) / det;
2175 inverse_jacobian(1, 2) =
2176 -(jacobian(0, 0) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 0)) /
2178 inverse_jacobian(2, 0) =
2179 (jacobian(1, 0) * jacobian(2, 1) - jacobian(1, 1) * jacobian(2, 0)) / det;
2180 inverse_jacobian(2, 1) =
2181 -(jacobian(0, 0) * jacobian(2, 1) - jacobian(0, 1) * jacobian(2, 0)) /
2183 inverse_jacobian(2, 2) =
2184 (jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0)) / det;
2201 const unsigned el_dim =
dim();
2207 return invert_jacobian<0>(jacobian, inverse_jacobian);
2210 return invert_jacobian<1>(jacobian, inverse_jacobian);
2213 return invert_jacobian<2>(jacobian, inverse_jacobian);
2216 return invert_jacobian<3>(jacobian, inverse_jacobian);
2220 std::ostringstream error_stream;
2221 error_stream <<
"Dimension of the element must be 0,1,2 or 3, not "
2222 << el_dim << std::endl;
2225 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2236 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<0>(
2242 oomph_info <<
"\nWarning: You are trying to calculate derivatives of "
2243 <<
"a jacobian w.r.t. nodal coordinates for a 'point' "
2244 <<
"element." << std::endl
2245 <<
"This makes no sense and is almost certainly an error."
2255 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<1>(
2261 const unsigned n_node =
nnode();
2264 for (
unsigned j = 0; j < n_node; j++)
2266 djacobian_dX(0, j) = dpsids(j, 0);
2275 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<2>(
2281 const unsigned n_node =
nnode();
2284 for (
unsigned j = 0; j < n_node; j++)
2287 djacobian_dX(0, j) =
2288 dpsids(j, 0) * jacobian(1, 1) - dpsids(j, 1) * jacobian(0, 1);
2291 djacobian_dX(1, j) =
2292 dpsids(j, 1) * jacobian(0, 0) - dpsids(j, 0) * jacobian(1, 0);
2301 void FiniteElement::dJ_eulerian_dnodal_coordinates_templated_helper<3>(
2307 const unsigned n_node =
nnode();
2310 for (
unsigned j = 0; j < n_node; j++)
2313 djacobian_dX(0, j) =
2315 (jacobian(1, 1) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 1)) +
2317 (jacobian(0, 2) * jacobian(2, 1) - jacobian(0, 1) * jacobian(2, 2)) +
2319 (jacobian(0, 1) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 1));
2322 djacobian_dX(1, j) =
2324 (jacobian(1, 2) * jacobian(2, 0) - jacobian(1, 0) * jacobian(2, 2)) +
2326 (jacobian(0, 0) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 0)) +
2328 (jacobian(0, 2) * jacobian(1, 0) - jacobian(0, 0) * jacobian(1, 2));
2331 djacobian_dX(2, j) =
2333 (jacobian(1, 0) * jacobian(2, 1) - jacobian(1, 1) * jacobian(2, 0)) +
2335 (jacobian(0, 1) * jacobian(2, 0) - jacobian(0, 0) * jacobian(2, 1)) +
2337 (jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0));
2347 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<0>(
2348 const double& det_jacobian,
2356 oomph_info <<
"\nWarning: You are trying to calculate derivatives of "
2357 <<
"eulerian derivatives of shape functions w.r.t. nodal "
2358 <<
"coordinates for a 'point' element." << std::endl
2359 <<
"This makes no sense and is almost certainly an error."
2370 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<1>(
2371 const double& det_jacobian,
2379 const double inv_det_jac = 1.0 / det_jacobian;
2382 const unsigned n_node =
nnode();
2385 for (
unsigned q = 0; q < n_node; q++)
2388 for (
unsigned j = 0; j < n_node; j++)
2390 d_dpsidx_dX(0, q, j, 0) =
2391 -djacobian_dX(0, q) * dpsids(j, 0) * inv_det_jac * inv_det_jac;
2402 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<2>(
2403 const double& det_jacobian,
2411 const double inv_det_jac = 1.0 / det_jacobian;
2414 const unsigned n_node =
nnode();
2417 for (
unsigned p = 0; p < 2; p++)
2420 for (
unsigned q = 0; q < n_node; q++)
2423 for (
unsigned j = 0; j < n_node; j++)
2426 d_dpsidx_dX(p, q, j, 0) =
2427 -djacobian_dX(p, q) * (inverse_jacobian(0, 0) * dpsids(j, 0) +
2428 inverse_jacobian(0, 1) * dpsids(j, 1));
2432 d_dpsidx_dX(p, q, j, 0) +=
2433 dpsids(j, 0) * dpsids(q, 1) - dpsids(j, 1) * dpsids(q, 0);
2435 d_dpsidx_dX(p, q, j, 0) *= inv_det_jac;
2438 d_dpsidx_dX(p, q, j, 1) =
2439 -djacobian_dX(p, q) * (inverse_jacobian(1, 1) * dpsids(j, 1) +
2440 inverse_jacobian(1, 0) * dpsids(j, 0));
2444 d_dpsidx_dX(p, q, j, 1) +=
2445 dpsids(j, 1) * dpsids(q, 0) - dpsids(j, 0) * dpsids(q, 1);
2447 d_dpsidx_dX(p, q, j, 1) *= inv_det_jac;
2459 void FiniteElement::d_dshape_eulerian_dnodal_coordinates_templated_helper<3>(
2460 const double& det_jacobian,
2468 const double inv_det_jac = 1.0 / det_jacobian;
2471 const unsigned n_node =
nnode();
2474 for (
unsigned p = 0; p < 3; p++)
2477 for (
unsigned q = 0; q < n_node; q++)
2480 for (
unsigned j = 0; j < n_node; j++)
2483 for (
unsigned i = 0;
i < 3;
i++)
2485 d_dpsidx_dX(p, q, j,
i) =
2486 -djacobian_dX(p, q) * (inverse_jacobian(
i, 0) * dpsids(j, 0) +
2487 inverse_jacobian(
i, 1) * dpsids(j, 1) +
2488 inverse_jacobian(
i, 2) * dpsids(j, 2));
2495 d_dpsidx_dX(p, q, j, 1) += ((dpsids(q, 2) * jacobian(1, 2) -
2496 dpsids(q, 1) * jacobian(2, 2)) *
2498 (dpsids(q, 0) * jacobian(2, 2) -
2499 dpsids(q, 2) * jacobian(0, 2)) *
2501 (dpsids(q, 1) * jacobian(0, 2) -
2502 dpsids(q, 0) * jacobian(1, 2)) *
2505 d_dpsidx_dX(p, q, j, 2) += ((dpsids(q, 1) * jacobian(2, 1) -
2506 dpsids(q, 2) * jacobian(1, 1)) *
2508 (dpsids(q, 2) * jacobian(0, 1) -
2509 dpsids(q, 0) * jacobian(2, 1)) *
2511 (dpsids(q, 0) * jacobian(1, 1) -
2512 dpsids(q, 1) * jacobian(0, 1)) *
2518 d_dpsidx_dX(p, q, j, 0) += ((dpsids(q, 1) * jacobian(2, 2) -
2519 dpsids(q, 2) * jacobian(1, 2)) *
2521 (dpsids(q, 2) * jacobian(0, 2) -
2522 dpsids(q, 0) * jacobian(2, 2)) *
2524 (dpsids(q, 0) * jacobian(1, 2) -
2525 dpsids(q, 1) * jacobian(0, 2)) *
2528 d_dpsidx_dX(p, q, j, 2) += ((dpsids(q, 2) * jacobian(1, 0) -
2529 dpsids(q, 1) * jacobian(2, 0)) *
2531 (dpsids(q, 0) * jacobian(2, 0) -
2532 dpsids(q, 2) * jacobian(0, 0)) *
2534 (dpsids(q, 1) * jacobian(0, 0) -
2535 dpsids(q, 0) * jacobian(1, 0)) *
2541 d_dpsidx_dX(p, q, j, 0) += ((dpsids(q, 2) * jacobian(1, 1) -
2542 dpsids(q, 1) * jacobian(2, 1)) *
2544 (dpsids(q, 0) * jacobian(2, 1) -
2545 dpsids(q, 2) * jacobian(0, 1)) *
2547 (dpsids(q, 1) * jacobian(0, 1) -
2548 dpsids(q, 0) * jacobian(1, 1)) *
2551 d_dpsidx_dX(p, q, j, 1) += ((dpsids(q, 1) * jacobian(2, 0) -
2552 dpsids(q, 2) * jacobian(1, 0)) *
2554 (dpsids(q, 2) * jacobian(0, 0) -
2555 dpsids(q, 0) * jacobian(2, 0)) *
2557 (dpsids(q, 0) * jacobian(1, 0) -
2558 dpsids(q, 1) * jacobian(0, 0)) *
2564 for (
unsigned i = 0;
i < 3;
i++)
2566 d_dpsidx_dX(p, q, j,
i) *= inv_det_jac;
2624 const unsigned el_dim =
dim();
2627 const unsigned n_shape =
nnode();
2634 std::ostringstream error_message;
2635 error_message <<
"Dimension mismatch" << std::endl;
2638 <<
" for the jacobian of the mapping to be well-defined"
2641 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2649 for (
unsigned i = 0;
i < el_dim;
i++)
2652 for (
unsigned j = 0; j < el_dim; j++)
2654 jacobian(
i, j) = 0.0;
2655 inverse_jacobian(
i, j) = 0.0;
2659 for (
unsigned l = 0; l < n_shape; l++)
2662 for (
unsigned k = 0; k < n_shape_type; k++)
2672 for (
unsigned i = 0;
i < el_dim;
i++)
2674 det *= jacobian(
i,
i);
2683 for (
unsigned i = 0;
i < el_dim;
i++)
2685 inverse_jacobian(
i,
i) = 1.0 / jacobian(
i,
i);
2705 const unsigned el_dim =
dim();
2709 const unsigned n_node =
nnode();
2712 if (djacobian_dX.
nrow() != el_dim)
2714 std::ostringstream error_message;
2715 error_message <<
"djacobian_dX must have the same number of rows as the"
2716 <<
"\nspatial dimension of the element.";
2718 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2721 if (djacobian_dX.
ncol() != n_node)
2723 std::ostringstream error_message;
2725 <<
"djacobian_dX must have the same number of columns as the"
2726 <<
"\nnumber of nodes in the element.";
2728 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2737 dJ_eulerian_dnodal_coordinates_templated_helper<0>(
2738 jacobian, dpsids, djacobian_dX);
2741 dJ_eulerian_dnodal_coordinates_templated_helper<1>(
2742 jacobian, dpsids, djacobian_dX);
2745 dJ_eulerian_dnodal_coordinates_templated_helper<2>(
2746 jacobian, dpsids, djacobian_dX);
2749 dJ_eulerian_dnodal_coordinates_templated_helper<3>(
2750 jacobian, dpsids, djacobian_dX);
2754 std::ostringstream error_stream;
2755 error_stream <<
"Dimension of the element must be 0,1,2 or 3, not "
2756 << el_dim << std::endl;
2759 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2780 const double& det_jacobian,
2788 const unsigned el_dim =
dim();
2792 const unsigned n_node =
nnode();
2795 if (d_dpsidx_dX.
nindex1() != el_dim || d_dpsidx_dX.
nindex2() != n_node ||
2796 d_dpsidx_dX.
nindex3() != n_node || d_dpsidx_dX.
nindex4() != el_dim)
2798 std::ostringstream error_message;
2799 error_message <<
"d_dpsidx_dX must be of the following dimensions:"
2800 <<
"\nd_dpsidx_dX(el_dim,n_node,n_node,el_dim)";
2802 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2811 d_dshape_eulerian_dnodal_coordinates_templated_helper<0>(
2820 d_dshape_eulerian_dnodal_coordinates_templated_helper<1>(
2829 d_dshape_eulerian_dnodal_coordinates_templated_helper<2>(
2838 d_dshape_eulerian_dnodal_coordinates_templated_helper<3>(
2848 std::ostringstream error_stream;
2849 error_stream <<
"Dimension of the element must be 0,1,2 or 3, not "
2850 << el_dim << std::endl;
2853 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2867 const unsigned n_basis = dbasis.
nindex1();
2868 const unsigned n_basis_type = dbasis.
nindex2();
2870 const unsigned n_dim =
dim();
2873 double new_derivatives[n_dim];
2876 for (
unsigned l = 0; l < n_basis; l++)
2879 for (
unsigned k = 0; k < n_basis_type; k++)
2882 for (
unsigned j = 0; j < n_dim; j++)
2885 new_derivatives[j] = 0.0;
2887 for (
unsigned i = 0;
i < n_dim;
i++)
2889 new_derivatives[j] += inverse_jacobian(j,
i) * dbasis(l, k,
i);
2893 for (
unsigned j = 0; j < n_dim; j++)
2895 dbasis(l, k, j) = new_derivatives[j];
2911 const unsigned n_basis = dbasis.
nindex1();
2912 const unsigned n_basis_type = dbasis.
nindex2();
2914 const unsigned n_dim =
dim();
2917 for (
unsigned l = 0; l < n_basis; l++)
2920 for (
unsigned k = 0; k < n_basis_type; k++)
2923 for (
unsigned j = 0; j < n_dim; j++)
2925 dbasis(l, k, j) *= inverse_jacobian(j, j);
2939 void FiniteElement::transform_second_derivatives_template<1>(
2947 const unsigned n_basis = dbasis.
nindex1();
2948 const unsigned n_basis_type = dbasis.
nindex2();
2952 for (
unsigned l = 0; l < n_basis; l++)
2955 for (
unsigned k = 0; k < n_basis_type; k++)
2958 d2basis(l, k, 0) / (jacobian(0, 0) * jacobian(0, 0))
2960 - dbasis(l, k, 0) * jacobian2(0, 0) /
2961 (jacobian(0, 0) * jacobian(0, 0) * jacobian(0, 0));
2978 void FiniteElement::transform_second_derivatives_template<2>(
2986 const unsigned n_basis = dbasis.
nindex1();
2987 const unsigned n_basis_type = dbasis.
nindex2();
2991 jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0);
3000 jacobian2(0, 0) * jacobian(1, 1) + jacobian(0, 0) * jacobian2(2, 1) -
3001 jacobian2(0, 1) * jacobian(1, 0) - jacobian(0, 1) * jacobian2(2, 0);
3003 jacobian2(2, 0) * jacobian(1, 1) + jacobian(0, 0) * jacobian2(1, 1) -
3004 jacobian2(2, 1) * jacobian(1, 0) - jacobian(0, 1) * jacobian2(1, 0);
3008 double dinverse_jacobiands[2][2][2];
3010 dinverse_jacobiands[0][0][0] =
3011 jacobian2(2, 1) / det - jacobian(1, 1) * ddetds[0] / (det * det);
3012 dinverse_jacobiands[0][1][0] =
3013 -jacobian2(0, 1) / det + jacobian(0, 1) * ddetds[0] / (det * det);
3014 dinverse_jacobiands[1][0][0] =
3015 -jacobian2(2, 0) / det + jacobian(1, 0) * ddetds[0] / (det * det);
3016 dinverse_jacobiands[1][1][0] =
3017 jacobian2(0, 0) / det - jacobian(0, 0) * ddetds[0] / (det * det);
3019 dinverse_jacobiands[0][0][1] =
3020 jacobian2(1, 1) / det - jacobian(1, 1) * ddetds[1] / (det * det);
3021 dinverse_jacobiands[0][1][1] =
3022 -jacobian2(2, 1) / det + jacobian(0, 1) * ddetds[1] / (det * det);
3023 dinverse_jacobiands[1][0][1] =
3024 -jacobian2(1, 0) / det + jacobian(1, 0) * ddetds[1] / (det * det);
3025 dinverse_jacobiands[1][1][1] =
3026 jacobian2(2, 0) / det - jacobian(0, 0) * ddetds[1] / (det * det);
3029 DShape dphidXds0(n_basis, n_basis_type, 2),
3030 dphidXds1(n_basis, n_basis_type, 2);
3032 for (
unsigned l = 0; l < n_basis; l++)
3034 for (
unsigned k = 0; k < n_basis_type; k++)
3036 for (
unsigned j = 0; j < 2; j++)
3040 dphidXds0(l, k, j) = dinverse_jacobiands[j][0][0] * dbasis(l, k, 0) +
3041 dinverse_jacobiands[j][1][0] * dbasis(l, k, 1) +
3042 inverse_jacobian(j, 0) * d2basis(l, k, 0) +
3043 inverse_jacobian(j, 1) * d2basis(l, k, 2);
3045 dphidXds1(l, k, j) = dinverse_jacobiands[j][0][1] * dbasis(l, k, 0) +
3046 dinverse_jacobiands[j][1][1] * dbasis(l, k, 1) +
3047 inverse_jacobian(j, 0) * d2basis(l, k, 2) +
3048 inverse_jacobian(j, 1) * d2basis(l, k, 1);
3054 for (
unsigned l = 0; l < n_basis; l++)
3056 for (
unsigned k = 0; k < n_basis_type; k++)
3059 for (
unsigned j = 0; j < 3; j++)
3061 d2basis(l, k, j) = 0.0;
3065 for (
unsigned i = 0;
i < 2;
i++)
3067 d2basis(l, k,
i) = inverse_jacobian(
i, 0) * dphidXds0(l, k,
i) +
3068 inverse_jacobian(
i, 1) * dphidXds1(l, k,
i);
3071 d2basis(l, k, 2) += inverse_jacobian(0, 0) * dphidXds0(l, k, 1) +
3072 inverse_jacobian(0, 1) * dphidXds1(l, k, 1);
3090 void FiniteElement::transform_second_derivatives_diagonal<1>(
3097 FiniteElement::transform_second_derivatives_template<1>(
3098 jacobian, inverse_jacobian, jacobian2, dbasis, d2basis);
3108 void FiniteElement::transform_second_derivatives_diagonal<2>(
3116 const unsigned n_basis = dbasis.
nindex1();
3117 const unsigned n_basis_type = dbasis.
nindex2();
3123 for (
unsigned l = 0; l < n_basis; l++)
3125 for (
unsigned k = 0; k < n_basis_type; k++)
3129 d2basis(l, k, 0) / (jacobian(0, 0) * jacobian(0, 0)) -
3130 dbasis(l, k, 0) * jacobian2(0, 0) /
3131 (jacobian(0, 0) * jacobian(0, 0) * jacobian(0, 0));
3134 d2basis(l, k, 1) / (jacobian(1, 1) * jacobian(1, 1)) -
3135 dbasis(l, k, 1) * jacobian2(1, 1) /
3136 (jacobian(1, 1) * jacobian(1, 1) * jacobian(1, 1));
3138 d2basis(l, k, 2) = d2basis(l, k, 2) / (jacobian(0, 0) * jacobian(1, 1));
3163 const unsigned el_dim =
dim();
3168 transform_second_derivatives_template<1>(
3169 jacobian, inverse_jacobian, jacobian2, dbasis, d2basis);
3172 transform_second_derivatives_template<2>(
3173 jacobian, inverse_jacobian, jacobian2, dbasis, d2basis);
3177 throw OomphLibError(
"Not implemented yet ... maybe one day",
3178 OOMPH_CURRENT_FUNCTION,
3179 OOMPH_EXCEPTION_LOCATION);
3187 std::ostringstream error_stream;
3188 error_stream <<
"Dimension of the element must be 1,2 or 3, not "
3189 << el_dim << std::endl;
3192 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3227 unsigned n_coordinates = s_fraction.size();
3228 for (
unsigned i = 0;
i < n_coordinates;
i++)
3253 const unsigned el_dim =
dim();
3257 for (
unsigned i = 0;
i < el_dim;
i++)
3274 const unsigned el_dim =
dim();
3278 for (
unsigned i = 0;
i < el_dim;
i++)
3310 const unsigned el_dim =
dim();
3314 for (
unsigned i = 0;
i < el_dim;
i++)
3333 const unsigned el_dim =
dim();
3360 const unsigned el_dim =
dim();
3388 const unsigned el_dim =
dim();
3391 unsigned nnod =
nnode();
3392 DShape dpsi(nnod, el_dim);
3421 const unsigned& ipt,
3428 const unsigned el_dim =
dim();
3449 det, jacobian, djacobian_dX, inverse_jacobian, dpsi, d_dpsidx_dX);
3486 const unsigned el_dim =
dim();
3488 const unsigned n_deriv =
N2deriv[el_dim];
3506 jacobian, inverse_jacobian, jacobian2, dpsi, d2psi);
3540 const unsigned el_dim =
dim();
3542 const unsigned n_deriv =
N2deriv[el_dim];
3560 jacobian, inverse_jacobian, jacobian2, dpsi, d2psi);
3575 const bool& store_local_dof_pt)
3578 const unsigned n_node =
nnode();
3589 for (
unsigned n = 1; n < n_node; n++)
3603 if (n_total_values == 0)
3615 for (
unsigned i = 0;
i < n_total_values;
i++)
3621 for (
unsigned n = 1; n < n_node; ++n)
3633 std::deque<unsigned long> global_eqn_number_queue;
3636 for (
unsigned n = 0; n < n_node; n++)
3642 unsigned n_value = nod_pt->
nvalue();
3645 for (
unsigned j = 0; j < n_value; j++)
3653 global_eqn_number_queue.push_back(
eqn_number);
3655 if (store_local_dof_pt)
3676 if (store_local_dof_pt)
3694 const unsigned n_node =
nnode();
3706 const unsigned n_dof =
ndof();
3711 int local_unknown = 0;
3717 for (
unsigned n = 0; n < n_node; n++)
3723 for (
unsigned i = 0;
i < n_value;
i++)
3728 if (local_unknown >= 0)
3734 const double old_var = *value_pt;
3737 *value_pt += fd_step;
3746 for (
unsigned m = 0; m < n_dof; m++)
3748 double sum = (newres[m] - residuals[m]) / fd_step;
3750 jacobian(m, local_unknown) = sum;
3754 *value_pt = old_var;
3778 unsigned n_nod =
nnode();
3781 if (n_nod == 0)
return;
3787 unsigned n_dof =
ndof();
3798 for (
unsigned j = 0; j < n_nod; j++)
3804 for (
unsigned i = 0;
i < dim_nod;
i++)
3807 double backup = nod_pt->
x(
i);
3811 nod_pt->
x(
i) += eps_fd;
3823 for (
unsigned l = 0; l < n_dof; l++)
3825 dresidual_dnodal_coordinates(l,
i, j) =
3826 (res_pls[l] - res[l]) / eps_fd;
3831 nod_pt->
x(
i) = backup;
3849 unsigned n_node =
nnode();
3855 std::vector<int> local_node_number;
3857 for (
unsigned i = 0;
i < n_node;
i++)
3864 local_node_number.push_back(
i);
3871 std::ostringstream error_stream;
3872 error_stream <<
"Node " << global_node_pt <<
" appears " << count
3873 <<
" times in an element." << std::endl
3874 <<
"In positions: ";
3875 for (std::vector<int>::iterator it = local_node_number.begin();
3876 it != local_node_number.end();
3879 error_stream << *it <<
" ";
3881 error_stream << std::endl <<
"That seems very odd." << std::endl;
3884 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
3890 for (
unsigned i = 0;
i < n_node;
i++)
3920 const unsigned n_node =
Nnode;
3922 for (
unsigned n = 0; n < n_node; n++)
3927 for (
unsigned i = 0;
i < el_dim;
i++)
3932 if (std::fabs(
s[
i] - s_node[
i]) > tol)
3962 unsigned nnod =
nnode();
3963 for (
unsigned j = 0; j < nnod; j++)
3972 cog[
i] /= double(nnod);
3976 for (
unsigned j = 0; j < nnod; j++)
3978 double dist_squared = 0.0;
3984 if (dist_squared > max_radius) max_radius = dist_squared;
3986 max_radius = sqrt(max_radius);
3993 const unsigned&
i)
const
3996 const unsigned n_node =
nnode();
4000 Shape psi(n_node, n_position_type);
4007 for (
unsigned l = 0; l < n_node; l++)
4010 for (
unsigned k = 0; k < n_position_type; k++)
4025 const unsigned&
i)
const
4028 const unsigned n_node =
nnode();
4033 Shape psi(n_node, n_position_type);
4040 for (
unsigned l = 0; l < n_node; l++)
4043 for (
unsigned k = 0; k < n_position_type; k++)
4059 const unsigned n_node =
nnode();
4066 Shape psi(n_node, n_position_type);
4071 for (
unsigned i = 0;
i < nodal_dim;
i++)
4076 for (
unsigned l = 0; l < n_node; l++)
4079 for (
unsigned k = 0; k < n_position_type; k++)
4096 const unsigned n_node =
nnode();
4103 Shape psi(n_node, n_position_type);
4108 for (
unsigned i = 0;
i < nodal_dim;
i++)
4113 for (
unsigned l = 0; l < n_node; l++)
4116 for (
unsigned k = 0; k < n_position_type; k++)
4136 const unsigned n_node =
nnode();
4140 const unsigned n_dim_element =
dim();
4143 Shape psi(n_node, n_position_type);
4144 DShape dpsids(n_node, n_position_type, n_dim_element);
4154 for (
unsigned i = 0;
i < n_dim_element;
i++)
4156 for (
unsigned j = 0; j < n_dim_element; j++)
4158 for (
unsigned k = 0; k < n_dim_node; k++)
4160 G(
i, j) += interpolated_G(
i, k) * interpolated_G(j, k);
4167 switch (n_dim_element)
4170 throw OomphLibError(
"Cannot calculate J_eulerian() for point element\n",
4171 OOMPH_CURRENT_FUNCTION,
4172 OOMPH_EXCEPTION_LOCATION);
4178 det = G(0, 0) * G(1, 1) - G(0, 1) * G(1, 0);
4181 det = G(0, 0) * G(1, 1) * G(2, 2) + G(0, 1) * G(1, 2) * G(2, 0) +
4182 G(0, 2) * G(1, 0) * G(2, 1) - G(0, 0) * G(1, 2) * G(2, 1) -
4183 G(0, 1) * G(1, 0) * G(2, 2) - G(0, 2) * G(1, 1) * G(2, 0);
4186 oomph_info <<
"More than 3 dimensions in J_eulerian()" << std::endl;
4201 const unsigned n_node =
nnode();
4206 const unsigned n_dim_element =
dim();
4209 Shape psi(n_node, n_position_type);
4210 DShape dpsids(n_node, n_position_type, n_dim_element);
4223 for (
unsigned i = 0;
i < n_dim_element;
i++)
4225 for (
unsigned j = 0; j < n_dim_element; j++)
4227 for (
unsigned k = 0; k < n_dim_node; k++)
4229 G(
i, j) += interpolated_G(
i, k) * interpolated_G(j, k);
4236 switch (n_dim_element)
4239 throw OomphLibError(
"Cannot calculate J_eulerian() for point element\n",
4240 OOMPH_CURRENT_FUNCTION,
4241 OOMPH_EXCEPTION_LOCATION);
4247 det = G(0, 0) * G(1, 1) - G(0, 1) * G(1, 0);
4250 det = G(0, 0) * G(1, 1) * G(2, 2) + G(0, 1) * G(1, 2) * G(2, 0) +
4251 G(0, 2) * G(1, 0) * G(2, 1) - G(0, 0) * G(1, 2) * G(2, 1) -
4252 G(0, 1) * G(1, 0) * G(2, 2) - G(0, 2) * G(1, 1) * G(2, 0);
4255 oomph_info <<
"More than 3 dimensions in J_eulerian()" << std::endl;
4277 const unsigned n_node =
nnode();
4286 const unsigned n_dim_element =
dim();
4289 Shape psi(n_node, n_position_type);
4290 DShape dpsi(n_node, n_dim_element);
4293 for (
unsigned ipt = 0; ipt < nintpt; ipt++)
4327 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
4351 const unsigned ncomponents = integral.size();
4353 for (
unsigned i = 0;
i < ncomponents;
i++)
4369 const unsigned n_dim =
dim();
4373 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
4376 for (
unsigned i = 0;
i < n_dim;
i++)
4382 for (
unsigned i = 0;
i < n_dim_eulerian;
i++)
4394 integrand_fct_pt(time, x, integrand);
4397 for (
unsigned i = 0;
i < ncomponents;
i++)
4399 integral[
i] += integrand[
i] * w * J;
4412 const unsigned ncomponents = integral.size();
4414 for (
unsigned i = 0;
i < ncomponents;
i++)
4430 const unsigned n_dim =
dim();
4434 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
4437 for (
unsigned i = 0;
i < n_dim;
i++)
4443 for (
unsigned i = 0;
i < n_dim_eulerian;
i++)
4455 integrand_fct_pt(x, integrand);
4458 for (
unsigned i = 0;
i < ncomponents;
i++)
4460 integral[
i] += integrand[
i] * w * J;
4485 OomphLibWarning(
"Pointer to spatial integration scheme has not been set.",
4486 "FiniteElement::self_test()",
4487 OOMPH_EXCEPTION_LOCATION);
4492 const unsigned dim_el =
dim();
4506 const unsigned n_node =
nnode();
4509 Shape psi(n_node, n_position_type);
4510 DShape dpsidx(n_node, dim_el);
4523 if (face_el_pt == 0)
4533 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
4536 for (
unsigned i = 0;
i < dim_el;
i++)
4543 for (
unsigned test = 0; test < ntest; test++)
4567 OOMPH_CURRENT_FUNCTION,
4568 OOMPH_EXCEPTION_LOCATION);
4573 if (std::fabs(jacobian) < 1.0e-16)
4575 std::ostringstream warning_stream;
4576 warning_stream <<
"Determinant of Jacobian matrix is zero at ipt "
4577 << ipt << std::endl;
4579 "FiniteElement::self_test()",
4580 OOMPH_EXCEPTION_LOCATION);
4589 std::ostringstream warning_stream;
4590 warning_stream <<
"Jacobian negative at integration point ipt="
4591 << ipt << std::endl;
4593 <<
"If you think that this is what you want you may: "
4595 <<
"set the (static) flag "
4596 <<
"FiniteElement::Accept_negative_jacobian to be true"
4600 "FiniteElement::self_test()",
4601 OOMPH_EXCEPTION_LOCATION);
4628 const unsigned& t_deriv)
4631 const unsigned n_node =
nnode();
4634 Shape psi(n_node, n_position_type);
4642 for (
unsigned l = 0; l < n_node; l++)
4645 for (
unsigned k = 0; k < n_position_type; k++)
4660 const unsigned& t_deriv,
4664 const unsigned n_node =
nnode();
4669 Shape psi(n_node, n_position_type);
4673 for (
unsigned i = 0;
i < nodal_dim;
i++)
4680 for (
unsigned l = 0; l < n_node; l++)
4683 for (
unsigned k = 0; k < n_position_type; k++)
4717 const unsigned n_node = this->
nnode();
4721 Shape psi(n_node, n_position_type);
4723 this->
shape(s, psi);
4726 const unsigned ncoord = this->
dim();
4728 for (
unsigned i = 0;
i < ncoord;
i++)
4735 for (
unsigned l = 0; l < n_node; l++)
4737 for (
unsigned k = 0; k < n_position_type; k++)
4740 const double psi_ = psi(l, k);
4741 for (
unsigned i = 0;
i < ncoord;
i++)
4767 const bool& use_coordinate_as_initial_guess)
4771 unsigned ncoord = this->
dim();
4779 double max_radius = 0.0;
4783 double radius = 0.0;
4784 for (
unsigned i = 0;
i < ncoord;
i++)
4786 radius += (cog[
i] - zeta[
i]) * (cog[
i] - zeta[
i]);
4788 radius = sqrt(radius);
4789 if (radius > Locate_zeta_helpers::
4812 if (!use_coordinate_as_initial_guess)
4818 unsigned num_plot_points =
4820 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
4824 s_list.push_back(s_c);
4832 bool keep_going =
true;
4838 if (use_coordinate_as_initial_guess)
4844 for (
unsigned i = 0;
i < ncoord;
i++)
4846 dx[
i] = zeta[
i] - inter_x[
i];
4852 double my_min_resid = DBL_MAX;
4857 unsigned n_list_coord = s_list.size();
4859 for (
unsigned i_coord = 0; i_coord < n_list_coord; i_coord++)
4861 for (
unsigned i = 0;
i < ncoord;
i++)
4863 s_local[
i] = s_list[i_coord][
i];
4869 for (
unsigned i = 0;
i < ncoord;
i++)
4871 work_dx[
i] = zeta[
i] - work_x[
i];
4874 double maxres = std::fabs(
4875 *std::max_element(work_dx.begin(), work_dx.end(),
AbsCmp<double>()));
4878 if (maxres < my_min_resid)
4880 my_min_resid = maxres;
4905 std::fabs(*std::max_element(dx.begin(), dx.end(),
AbsCmp<double>()));
4929 for (
unsigned i = 0;
i < ncoord;
i++)
4932 work_s[
i] += fd_step;
4938 for (
unsigned j = 0; j < ncoord; j++)
4940 jacobian(j,
i) = -(work_r[j] - r[j]) / fd_step;
4950 unsigned n_node = this->
nnode();
4952 Shape psi(n_node, n_position_type);
4953 DShape dpsids(n_node, n_position_type, ncoord);
4982 for (
unsigned l = 0; l < n_node; l++)
4986 for (
unsigned k = 0; k < n_position_type; k++)
4989 for (
unsigned i = 0;
i < ncoord;
i++)
4991 for (
unsigned j = 0; j < ncoord; j++)
4993 interpolated_dxds(
i, j) +=
5002 for (
unsigned i = 0;
i < ncoord;
i++)
5004 for (
unsigned j = 0; j < ncoord; j++)
5006 jacobian(
i, j) = -interpolated_dxds(
i, j);
5022 <<
"FiniteElement::locate_zeta()" << std::endl;
5023 oomph_info <<
"Should not affect the result!" << std::endl;
5028 for (
unsigned i = 0;
i < ncoord;
i++)
5035 for (
unsigned i = 0;
i < ncoord;
i++)
5037 dx[
i] = zeta[
i] - inter_x[
i];
5042 std::fabs(*std::max_element(dx.begin(), dx.end(),
AbsCmp<double>()));
5050 }
while (keep_going);
5063 for (
unsigned i = 0;
i < ncoord;
i++)
5065 dx[
i] = zeta[
i] - inter_x[
i];
5070 std::fabs(*std::max_element(dx.begin(), dx.end(),
AbsCmp<double>()));
5094 geom_object_pt =
this;
5104 const unsigned n_node =
nnode();
5105 for (
unsigned n = 0; n < n_node; n++)
5127 std::set<std::pair<Data*, unsigned>>& paired_field_data)
5131 for (
unsigned n = 0; n < n_internal; n++)
5136 const unsigned n_value = dat_pt->
nvalue();
5139 for (
unsigned i = 0;
i < n_value;
i++)
5141 paired_field_data.insert(std::make_pair(dat_pt,
i));
5146 const unsigned n_node = this->
nnode();
5147 for (
unsigned n = 0; n < n_node; n++)
5152 const unsigned n_value = nod_pt->
nvalue();
5155 for (
unsigned i = 0;
i < n_value;
i++)
5157 paired_field_data.insert(std::make_pair(nod_pt,
i));
5171 #ifdef OOMPH_HAS_MPI
5201 for (
unsigned i = 0;
i < nnode_face;
i++)
5233 unsigned n_dim =
dim();
5241 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
5247 for (
unsigned i = 0;
i < n_dim + 1;
i++)
5255 for (
unsigned i = 0;
i < n_dim;
i++)
5257 outfile << zeta[
i] <<
" ";
5259 outfile << std::endl;
5275 unsigned n_dim_el = this->
dim();
5279 if (n_dim_el == 0)
return 1.0;
5282 unsigned n_node =
nnode();
5291 Shape psi(n_node, n_position_type);
5292 DShape dpsids(n_node, n_position_type, n_dim_el);
5301 for (
unsigned l = 0; l < n_node; l++)
5304 for (
unsigned k = 0; k < n_position_type; k++)
5307 for (
unsigned i = 0;
i < n_dim;
i++)
5310 for (
unsigned j = 0; j < n_dim_el; j++)
5312 interpolated_A(j,
i) +=
5320 for (
unsigned i = 0;
i < n_dim_el;
i++)
5322 for (
unsigned j = 0; j < n_dim_el; j++)
5325 for (
unsigned k = 0; k < n_dim; k++)
5327 A(
i, j) += interpolated_A(
i, k) * interpolated_A(j, k);
5340 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
5344 "FaceElement::J_eulerian()",
5345 OOMPH_EXCEPTION_LOCATION);
5361 const unsigned n_node =
nnode();
5368 const unsigned n_dim_el =
dim();
5371 Shape psi(n_node, n_position_type);
5372 DShape dpsids(n_node, n_position_type, n_dim_el);
5384 for (
unsigned l = 0; l < n_node; l++)
5387 for (
unsigned k = 0; k < n_position_type; k++)
5390 for (
unsigned i = 0;
i < n_dim;
i++)
5393 for (
unsigned j = 0; j < n_dim_el; j++)
5395 interpolated_A(j,
i) +=
5404 for (
unsigned i = 0;
i < n_dim_el;
i++)
5406 for (
unsigned j = 0; j < n_dim_el; j++)
5409 for (
unsigned k = 0; k < n_dim; k++)
5411 A(
i, j) += interpolated_A(
i, k) * interpolated_A(j, k);
5424 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
5428 "FaceElement::J_eulerian_at_knot()",
5429 OOMPH_EXCEPTION_LOCATION);
5446 const unsigned n_node =
nnode();
5453 const unsigned n_dim_el =
dim();
5456 Shape psi(n_node, n_position_type);
5457 DShape dpsids(n_node, n_position_type, n_dim_el);
5461 for (
unsigned ipt = 0; ipt < nintpt; ipt++)
5473 for (
unsigned l = 0; l < n_node; l++)
5476 for (
unsigned k = 0; k < n_position_type; k++)
5479 for (
unsigned i = 0;
i < n_dim;
i++)
5483 for (
unsigned j = 0; j < n_dim_el; j++)
5485 interpolated_A(j,
i) +=
5495 for (
unsigned i = 0;
i < n_dim_el;
i++)
5497 for (
unsigned j = 0; j < n_dim_el; j++)
5500 for (
unsigned k = 0; k < n_dim; k++)
5502 A(
i, j) += interpolated_A(
i, k) * interpolated_A(j, k);
5515 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
5519 "FaceElement::J_eulerian_at_knot()",
5520 OOMPH_EXCEPTION_LOCATION);
5548 const unsigned element_dim =
dim();
5556 if (
s.size() != element_dim)
5558 std::ostringstream error_stream;
5560 <<
"Local coordinate s passed to outer_unit_normal() has dimension "
5561 <<
s.size() << std::endl
5562 <<
"but element dimension is " << element_dim << std::endl;
5565 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5573 std::ostringstream error_stream;
5574 error_stream <<
"Tangent direction vector has dimension "
5576 <<
"but spatial dimension is " << spatial_dim << std::endl;
5579 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5583 if (unit_normal.size() != spatial_dim)
5585 std::ostringstream error_stream;
5586 error_stream <<
"Unit normal passed to outer_unit_normal() has dimension "
5587 << unit_normal.size() << std::endl
5588 <<
"but spatial dimension is " << spatial_dim << std::endl;
5591 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5596 unsigned ntangent_vec = tang_vec.size();
5602 switch (element_dim)
5608 if (ntangent_vec != 1)
5610 std::ostringstream error_stream;
5612 <<
"This is a 0D FaceElement, we need one tangent vector.\n"
5613 <<
"You have given me " << tang_vec.size() <<
" vector(s).\n";
5615 OOMPH_CURRENT_FUNCTION,
5616 OOMPH_EXCEPTION_LOCATION);
5625 if (ntangent_vec != 1)
5627 std::ostringstream error_stream;
5629 <<
"This is a 1D FaceElement, we need one tangent vector.\n"
5630 <<
"You have given me " << tang_vec.size() <<
" vector(s).\n";
5632 OOMPH_CURRENT_FUNCTION,
5633 OOMPH_EXCEPTION_LOCATION);
5642 if (ntangent_vec != 2)
5644 std::ostringstream error_stream;
5646 <<
"This is a 2D FaceElement, we need two tangent vectors.\n"
5647 <<
"You have given me " << tang_vec.size() <<
" vector(s).\n";
5649 OOMPH_CURRENT_FUNCTION,
5650 OOMPH_EXCEPTION_LOCATION);
5658 "Cannot have a FaceElement with dimension higher than 2",
5659 OOMPH_CURRENT_FUNCTION,
5660 OOMPH_EXCEPTION_LOCATION);
5665 for (
unsigned vec_i = 0; vec_i < ntangent_vec; vec_i++)
5667 if (tang_vec[vec_i].
size() != spatial_dim)
5669 std::ostringstream error_stream;
5670 error_stream <<
"This problem has " << spatial_dim
5671 <<
" spatial dimensions.\n"
5672 <<
"But the " << vec_i <<
" tangent vector has size "
5673 << tang_vec[vec_i].size() << std::endl;
5675 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5683 switch (element_dim)
5688 std::ostringstream error_stream;
5690 <<
"It is unclear how to calculate a tangent and normal vectors for "
5691 <<
"a point element.\n";
5693 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
5706 const unsigned n_position_type_bulk =
5716 Shape psi(n_node_bulk, n_position_type_bulk);
5717 DShape dpsids(n_node_bulk, n_position_type_bulk, 2);
5726 for (
unsigned l = 0; l < n_node_bulk; l++)
5729 for (
unsigned k = 0; k < n_position_type_bulk; k++)
5732 for (
unsigned j = 0; j < 2; j++)
5735 for (
unsigned i = 0;
i < spatial_dim;
i++)
5738 interpolated_dxds(j,
i) +=
5758 unsigned interior_direction = 0;
5761 for (
unsigned i = 0;
i < spatial_dim;
i++)
5765 tangent[
i] = interpolated_dxds(0,
i) * dsbulk_dsface(0, 0) +
5766 interpolated_dxds(1,
i) * dsbulk_dsface(1, 0);
5769 interior_tangent[
i] = interpolated_dxds(interior_direction,
i);
5783 for (
unsigned i = 0;
i < spatial_dim;
i++)
5785 unit_normal[
i] = full_normal[
i];
5790 for (
unsigned i = 0;
i < spatial_dim;
i++)
5796 tang_vec[0][0] = -unit_normal[1];
5797 tang_vec[0][1] = unit_normal[0];
5811 if (spatial_dim != 3)
5813 std::ostringstream error_stream;
5814 error_stream <<
"There are only " << spatial_dim
5815 <<
"coordinates at the nodes of this 2D FaceElement,\n"
5816 <<
"which must have come from a 3D Bulk element\n";
5818 OOMPH_CURRENT_FUNCTION,
5819 OOMPH_EXCEPTION_LOCATION);
5824 const unsigned n_node = this->
nnode();
5830 Shape psi(n_node, n_position_type);
5831 DShape dpsids(n_node, n_position_type, 2);
5840 for (
unsigned l = 0; l < n_node; l++)
5843 for (
unsigned k = 0; k < n_position_type; k++)
5846 for (
unsigned j = 0; j < 2; j++)
5849 for (
unsigned i = 0;
i < 3;
i++)
5854 interpolated_dxds[j][
i] +=
5864 interpolated_dxds[0], interpolated_dxds[1], unit_normal);
5869 for (
unsigned i = 0;
i < spatial_dim;
i++)
5894 std::ostringstream err_stream;
5895 err_stream <<
"The angle between the unit normal and the \n"
5896 <<
"direction vector is less than 10 degrees.";
5898 OOMPH_CURRENT_FUNCTION,
5899 OOMPH_EXCEPTION_LOCATION);
5913 tang_vec[0][0] = t1_scaling * interpolated_dxds[0][0] +
5914 t2_scaling * interpolated_dxds[1][0];
5915 tang_vec[0][1] = t1_scaling * interpolated_dxds[0][1] +
5916 t2_scaling * interpolated_dxds[1][1];
5917 tang_vec[0][2] = t1_scaling * interpolated_dxds[0][2] +
5918 t2_scaling * interpolated_dxds[1][2];
5925 for (
unsigned vec_i = 0; vec_i < 2; vec_i++)
5930 for (
unsigned dim_i = 0; dim_i < spatial_dim; dim_i++)
5932 tang_vec[vec_i][dim_i] /= tang_length;
5949 if (a != 0.0 || b != 0.0)
5951 double a_sq = a * a;
5952 double b_sq = b * b;
5953 double c_sq = c * c;
5955 tang_vec[0][0] = -b / sqrt(a_sq + b_sq);
5956 tang_vec[0][1] = a / sqrt(a_sq + b_sq);
5959 double z = (a_sq + b_sq) / sqrt(a_sq * c_sq + b_sq * c_sq +
5960 (a_sq + b_sq) * (a_sq + b_sq));
5962 tang_vec[1][0] = -(a * c * z) / (a_sq + b_sq);
5963 tang_vec[1][1] = -(b * c * z) / (a_sq + b_sq);
5972 std::ostringstream warn_stream;
5974 <<
"I have detected that your normal vector is [0,0,1].\n"
5975 <<
"Since this function compares against 0.0, you may\n"
5976 <<
"get discontinuous tangent vectors.";
5978 OOMPH_CURRENT_FUNCTION,
5979 OOMPH_EXCEPTION_LOCATION);
5982 tang_vec[0][0] = 1.0;
5983 tang_vec[0][1] = 0.0;
5984 tang_vec[0][2] = 0.0;
5986 tang_vec[1][0] = 0.0;
5987 tang_vec[1][1] = 1.0;
5988 tang_vec[1][2] = 0.0;
5993 OOMPH_CURRENT_FUNCTION,
5994 OOMPH_EXCEPTION_LOCATION);
6003 "Cannot have a FaceElement with dimension higher than 2",
6004 OOMPH_CURRENT_FUNCTION,
6005 OOMPH_EXCEPTION_LOCATION);
6017 const unsigned& ipt,
6022 const unsigned element_dim =
dim();
6025 for (
unsigned i = 0;
i < element_dim;
i++)
6040 const unsigned element_dim =
dim();
6048 if (
s.size() != element_dim)
6050 std::ostringstream error_stream;
6052 <<
"Local coordinate s passed to outer_unit_normal() has dimension "
6053 <<
s.size() << std::endl
6054 <<
"but element dimension is " << element_dim << std::endl;
6057 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
6061 if (unit_normal.size() != spatial_dim)
6063 std::ostringstream error_stream;
6064 error_stream <<
"Unit normal passed to outer_unit_normal() has dimension "
6065 << unit_normal.size() << std::endl
6066 <<
"but spatial dimension is " << spatial_dim << std::endl;
6069 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);