613 std::vector<std::vector<double>> sList0;
615 std::vector<unsigned> nList;
621 if (std::is_base_of<SolidTElement<3, 2>,
ELEMENT>::value)
632 nList = {0, 1, 2, 3};
649 nList = {0, 4, 6, 2, 1, 5, 7, 3};
655 std::vector<Vector<double>> sList;
656 for (
auto s0 : sList0)
671 Vector<double> coordinate;
675 std::vector<double> value;
677 std::vector<Data> data;
679 std::vector<Point> points;
680 points.reserve(nel * sList.size());
685 std::vector<unsigned long> connectivity;
686 unsigned long offset;
689 std::vector<Cell> cells;
693 for (
unsigned e = 0; e < nel; e++)
696 auto el_pt =
dynamic_cast<ELEMENT*
>(
699 std::vector<unsigned long> connectivity;
701 for (
const auto& s : sList)
705 auto n_pt =
dynamic_cast<SolidNode*
>(el_pt->node_pt(
n));
708 el_pt->interpolated_x(s, x);
710 Vector<double> dxdt(3);
711 el_pt->interpolated_dxdt(s, 1, dxdt);
713 Vector<double> xi(3);
714 el_pt->interpolated_xi(s, xi);
716 Vector<double> body_force(3);
717 auto bodyForceFct = el_pt->body_force_fct_pt();
718 if (bodyForceFct) bodyForceFct(time(), xi, body_force);
720 DenseMatrix<double> sigma(3,3);
721 el_pt->get_stress(s, sigma);
722 double pressure = (sigma(0,0)+sigma(1,1)+sigma(2,2))/3;
724 DenseMatrix<double> strain(3,3);
725 el_pt->get_strain(s, strain);
727 std::vector<double> dudt
728 {el_pt->dnodal_position_dt(
n, 0),
729 el_pt->dnodal_position_dt(
n, 1),
730 el_pt->dnodal_position_dt(
n, 2)};
732 std::vector<double> coupledForce(3,0.0);
733 double couplingWeight = 0.0;
735 if (c_el_pt !=
nullptr) {
736 double nnode = c_el_pt->nnode();
737 double dim = c_el_pt->dim();
739 c_el_pt->shape(s, psi);
740 for (
int n=0;
n<nnode; ++
n) {
741 for (
int d=0; d<dim; ++d) {
742 coupledForce[d] += c_el_pt->get_coupling_residual(
n, d) * psi(
n);
744 couplingWeight += c_el_pt->get_coupling_weight(
n) * psi(
n);
748 std::set<unsigned>* boundaries_pt;
749 n_pt->get_boundaries_pt(boundaries_pt);
750 double b = boundaries_pt ? *boundaries_pt->begin() : -1;
751 std::vector<double> pin {(
double) n_pt->position_is_pinned(0),
752 (
double) n_pt->position_is_pinned(1),
753 (
double) n_pt->position_is_pinned(2)};
756 {
"Velocity", {dxdt[0], dxdt[1], dxdt[2]}},
757 {
"Displacement", {x[0] - xi[0], x[1] - xi[1], x[2] - xi[2]}},
758 {
"BodyForce", {body_force[0], body_force[1], body_force[2]}},
761 {
"CoupledForce",coupledForce},
762 {
"CouplingWeight", {couplingWeight}},
775 connectivity.push_back(points.size() - 1);
778 cells.push_back({connectivity, points.size(), vtkFormat});
782 static unsigned count = 0;
786 std::string vtkFileName =
name_ +
"FEM_" + std::to_string(count++) +
".vtu";
789 std::ofstream vtk(vtkFileName);
792 vtk <<
"<?xml version=\"1.0\"?>\n"
794 "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n"
795 "<UnstructuredGrid>\n"
796 "<Piece NumberOfPoints=\""
798 <<
"\" NumberOfCells=\""
803 "<DataArray type=\"Float32\" Name=\"Position\" NumberOfComponents=\"3\" format=\"ascii\">\n";
804 for (
auto point : points)
806 vtk << point.coordinate[0] <<
" " << point.coordinate[1] <<
" " << point.coordinate[2] <<
"\n";
808 vtk <<
"</DataArray>\n"
810 if (not points.empty())
812 vtk <<
"<PointData Vectors=\"vector\">\n";
813 for (
int i = 0;
i < points.front().data.size(); ++
i)
815 auto data = points.front().data[
i];
816 vtk << R
"(<DataArray type="Float32" Name=")"
817 << points.front().data[i].name
818 << R"(" NumberOfComponents=")"
819 << points.front().data[i].value.size()
820 << R"(" format="ascii">)" << "\n";
821 for (
const Point& point : points)
823 for (
const auto& value : point.data[
i].value)
831 vtk <<
"</PointData>\n"
835 "<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
836 for (
const Cell& cell : cells)
838 for (
auto point : cell.connectivity)
843 vtk <<
"</DataArray>\n"
844 "<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n";
845 for (
const Cell& cell : cells)
847 vtk << cell.offset <<
" ";
849 vtk <<
"</DataArray>\n"
850 "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n";
851 for (
const Cell& cell : cells)
853 vtk << cell.type <<
" ";
860 "</UnstructuredGrid>\n"
#define OOMPH_MPI_PROCESSOR_ID
Definition: SolidProblem.h:62
double getOomphTime() const
get function for current time
Definition: SolidProblem.h:279
Wrapper class for solid elements to be coupled with discrete solid particles on the surfaces.
Definition: ScaleCoupledElement.h:37