Skip to content

Commit

Permalink
Clean up recording types
Browse files Browse the repository at this point in the history
  • Loading branch information
@TLCFEM
TLCFEM committed Oct 31, 2023
1 parent bddacbd commit ca27417
Show file tree
Hide file tree
Showing 38 changed files with 1,016 additions and 1,113 deletions.
2 changes: 1 addition & 1 deletion Constraint/ConstraintParser.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -504,7 +504,7 @@ int create_new_criterion(const shared_ptr<DomainBase>& domain, istringstream& co
return SUANPAN_SUCCESS;
}

domain->insert(make_shared<MaxHistory>(tag, step_tag, to_list(type), limit));
domain->insert(make_shared<MaxHistory>(tag, step_tag, to_token(type), limit));

return SUANPAN_SUCCESS;
}
Expand Down
36 changes: 18 additions & 18 deletions Domain/Node.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -537,24 +537,6 @@ std::vector<vec> Node::record(const OutputType L) const {
else if(L == OutputType::U) data.push_back(current_displacement);
else if(L == OutputType::V) data.push_back(current_velocity);
else if(L == OutputType::A) data.push_back(current_acceleration);
else if(L == OutputType::U1) data.emplace_back(vec{current_displacement.n_elem >= 1 ? current_displacement(0) : 0.});
else if(L == OutputType::U2) data.emplace_back(vec{current_displacement.n_elem >= 2 ? current_displacement(1) : 0.});
else if(L == OutputType::U3) data.emplace_back(vec{current_displacement.n_elem >= 3 ? current_displacement(2) : 0.});
else if(L == OutputType::U4 || L == OutputType::UR1) data.emplace_back(vec{current_displacement.n_elem >= 4 ? current_displacement(3) : 0.});
else if(L == OutputType::U5 || L == OutputType::UR2) data.emplace_back(vec{current_displacement.n_elem >= 5 ? current_displacement(4) : 0.});
else if(L == OutputType::U6 || L == OutputType::UR3) data.emplace_back(vec{current_displacement.n_elem >= 6 ? current_displacement(5) : 0.});
else if(L == OutputType::V1) data.emplace_back(vec{current_velocity.n_elem >= 1 ? current_velocity(0) : 0.});
else if(L == OutputType::V2) data.emplace_back(vec{current_velocity.n_elem >= 2 ? current_velocity(1) : 0.});
else if(L == OutputType::V3) data.emplace_back(vec{current_velocity.n_elem >= 3 ? current_velocity(2) : 0.});
else if(L == OutputType::V4 || L == OutputType::VR1) data.emplace_back(vec{current_velocity.n_elem >= 4 ? current_velocity(3) : 0.});
else if(L == OutputType::V5 || L == OutputType::VR2) data.emplace_back(vec{current_velocity.n_elem >= 5 ? current_velocity(4) : 0.});
else if(L == OutputType::V6 || L == OutputType::VR3) data.emplace_back(vec{current_velocity.n_elem >= 6 ? current_velocity(5) : 0.});
else if(L == OutputType::A1) data.emplace_back(vec{current_acceleration.n_elem >= 1 ? current_acceleration(0) : 0.});
else if(L == OutputType::A2) data.emplace_back(vec{current_acceleration.n_elem >= 2 ? current_acceleration(1) : 0.});
else if(L == OutputType::A3) data.emplace_back(vec{current_acceleration.n_elem >= 3 ? current_acceleration(2) : 0.});
else if(L == OutputType::A4 || L == OutputType::AR1) data.emplace_back(vec{current_acceleration.n_elem >= 4 ? current_acceleration(3) : 0.});
else if(L == OutputType::A5 || L == OutputType::AR2) data.emplace_back(vec{current_acceleration.n_elem >= 5 ? current_acceleration(4) : 0.});
else if(L == OutputType::A6 || L == OutputType::AR3) data.emplace_back(vec{current_acceleration.n_elem >= 6 ? current_acceleration(5) : 0.});
else if(L == OutputType::RF1) data.emplace_back(vec{current_resistance.n_elem >= 1 ? current_resistance(0) : 0.});
else if(L == OutputType::RF2) data.emplace_back(vec{current_resistance.n_elem >= 2 ? current_resistance(1) : 0.});
else if(L == OutputType::RF3) data.emplace_back(vec{current_resistance.n_elem >= 3 ? current_resistance(2) : 0.});
Expand All @@ -573,6 +555,24 @@ std::vector<vec> Node::record(const OutputType L) const {
else if(L == OutputType::IF4 || L == OutputType::IM1) data.emplace_back(vec{current_inertial_force.n_elem >= 4 ? current_inertial_force(3) : 0.});
else if(L == OutputType::IF5 || L == OutputType::IM2) data.emplace_back(vec{current_inertial_force.n_elem >= 5 ? current_inertial_force(4) : 0.});
else if(L == OutputType::IF6 || L == OutputType::IM3) data.emplace_back(vec{current_inertial_force.n_elem >= 6 ? current_inertial_force(5) : 0.});
else if(L == OutputType::U1) data.emplace_back(vec{current_displacement.n_elem >= 1 ? current_displacement(0) : 0.});
else if(L == OutputType::U2) data.emplace_back(vec{current_displacement.n_elem >= 2 ? current_displacement(1) : 0.});
else if(L == OutputType::U3) data.emplace_back(vec{current_displacement.n_elem >= 3 ? current_displacement(2) : 0.});
else if(L == OutputType::U4 || L == OutputType::UR1) data.emplace_back(vec{current_displacement.n_elem >= 4 ? current_displacement(3) : 0.});
else if(L == OutputType::U5 || L == OutputType::UR2) data.emplace_back(vec{current_displacement.n_elem >= 5 ? current_displacement(4) : 0.});
else if(L == OutputType::U6 || L == OutputType::UR3) data.emplace_back(vec{current_displacement.n_elem >= 6 ? current_displacement(5) : 0.});
else if(L == OutputType::V1) data.emplace_back(vec{current_velocity.n_elem >= 1 ? current_velocity(0) : 0.});
else if(L == OutputType::V2) data.emplace_back(vec{current_velocity.n_elem >= 2 ? current_velocity(1) : 0.});
else if(L == OutputType::V3) data.emplace_back(vec{current_velocity.n_elem >= 3 ? current_velocity(2) : 0.});
else if(L == OutputType::V4 || L == OutputType::VR1) data.emplace_back(vec{current_velocity.n_elem >= 4 ? current_velocity(3) : 0.});
else if(L == OutputType::V5 || L == OutputType::VR2) data.emplace_back(vec{current_velocity.n_elem >= 5 ? current_velocity(4) : 0.});
else if(L == OutputType::V6 || L == OutputType::VR3) data.emplace_back(vec{current_velocity.n_elem >= 6 ? current_velocity(5) : 0.});
else if(L == OutputType::A1) data.emplace_back(vec{current_acceleration.n_elem >= 1 ? current_acceleration(0) : 0.});
else if(L == OutputType::A2) data.emplace_back(vec{current_acceleration.n_elem >= 2 ? current_acceleration(1) : 0.});
else if(L == OutputType::A3) data.emplace_back(vec{current_acceleration.n_elem >= 3 ? current_acceleration(2) : 0.});
else if(L == OutputType::A4 || L == OutputType::AR1) data.emplace_back(vec{current_acceleration.n_elem >= 4 ? current_acceleration(3) : 0.});
else if(L == OutputType::A5 || L == OutputType::AR2) data.emplace_back(vec{current_acceleration.n_elem >= 5 ? current_acceleration(4) : 0.});
else if(L == OutputType::A6 || L == OutputType::AR3) data.emplace_back(vec{current_acceleration.n_elem >= 6 ? current_acceleration(5) : 0.});

return data;
}
Expand Down
30 changes: 9 additions & 21 deletions Element/ElementParser.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1257,7 +1257,7 @@ void new_dkt3(unique_ptr<Element>& return_obj, istringstream& command) {
return;
}

unsigned num_ip = 3;
auto num_ip = 3u;
if(!get_optional_input(command, num_ip)) {
suanpan_error("A valid number of integration points is required.\n");
return;
Expand Down Expand Up @@ -1291,7 +1291,7 @@ void new_dkt4(unique_ptr<Element>& return_obj, istringstream& command) {
return;
}

unsigned num_ip = 3;
auto num_ip = 3u;
if(!get_optional_input(command, num_ip)) {
suanpan_error("A valid number of integration points is required.\n");
return;
Expand Down Expand Up @@ -1325,7 +1325,7 @@ void new_dkts3(unique_ptr<Element>& return_obj, istringstream& command) {
return;
}

unsigned num_ip = 3;
auto num_ip = 3u;
if(!get_optional_input(command, num_ip)) {
suanpan_error("A valid number of integration points is required.\n");
return;
Expand Down Expand Up @@ -1462,13 +1462,13 @@ void new_f21(unique_ptr<Element>& return_obj, istringstream& command) {
return;
}

unsigned int_pt = 6;
auto int_pt = 6u;
if(!get_optional_input(command, int_pt)) {
suanpan_error("A valid number of integration points is required.\n");
return;
}

unsigned nonlinear = 0;
auto nonlinear = 0u;
if(command.eof())
suanpan_debug("Linear geometry assumed.\n");
else if(!get_input(command, nonlinear))
Expand Down Expand Up @@ -1910,16 +1910,8 @@ void new_pcpedc(unique_ptr<Element>& return_obj, istringstream& command, const u
}

double alpha, n, k;
if(!get_optional_input(command, alpha)) {
suanpan_error("A valid alpha is required.\n");
return;
}
if(!get_optional_input(command, n)) {
suanpan_error("A valid porosity is required.\n");
return;
}
if(!get_optional_input(command, k)) {
suanpan_error("A valid permeability is required.\n");
if(!get_input(command, alpha, n, k)) {
suanpan_error("A valid parameter is required.\n");
return;
}

Expand Down Expand Up @@ -1952,12 +1944,8 @@ void new_pcpeuc(unique_ptr<Element>& return_obj, istringstream& command, const u
}

double alpha, n;
if(!get_optional_input(command, alpha)) {
suanpan_error("A valid alpha is required.\n");
return;
}
if(!get_optional_input(command, n)) {
suanpan_error("A valid porosity is required.\n");
if(!get_input(command, alpha, n)) {
suanpan_error("A valid parameter is required.\n");
return;
}

Expand Down
87 changes: 1 addition & 86 deletions Element/Membrane/Drilling/GCMQ.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,77 +24,6 @@
#include <Toolbox/tensor.h>
#include <Toolbox/utility.h>

/**
* \brief create converter for resultant forces
* \param E edge label
* \param T element thickness
* \param EC element coordinates
* \param IP integration plan
* \param TRANS transformation matrix from parent to global
*/
GCMQ::ResultantConverter::ResultantConverter(const Edge E, const double T, const mat& EC, const IntegrationPlan& IP, const mat& TRANS)
: direction_cosine(2, 3) {
const auto &X1 = IP(0, 0), &X2 = IP(1, 0);

vec node_i, node_j, pt_a, pt_b;

switch(E) {
case Edge::A:
node_i = EC.row(0).t();
node_j = EC.row(1).t();
pt_a = TRANS * form_stress_mode(X1, -1.);
pt_b = TRANS * form_stress_mode(X2, -1.);
break;
case Edge::B:
node_i = EC.row(1).t();
node_j = EC.row(2).t();
pt_a = TRANS * form_stress_mode(1., X1);
pt_b = TRANS * form_stress_mode(1., X2);
break;
case Edge::C:
node_i = EC.row(2).t();
node_j = EC.row(3).t();
pt_a = TRANS * form_stress_mode(X1, 1.);
pt_b = TRANS * form_stress_mode(X2, 1.);
break;
case Edge::D:
node_i = EC.row(3).t();
node_j = EC.row(0).t();
pt_a = TRANS * form_stress_mode(-1., X1);
pt_b = TRANS * form_stress_mode(-1., X2);
break;
}

const vec incre = node_j - node_i;

const auto edge_length = norm(incre);

const auto angle = 2. * atan2(incre(1), incre(0)) - datum::pi;

const auto sin_angle = sin(angle);
const auto cos_angle = cos(angle);

// transformation to local reference frame
direction_cosine(0, 0) = .5 + .5 * cos_angle;
direction_cosine(0, 1) = .5 - .5 * cos_angle;
direction_cosine(0, 2) = sin_angle;
direction_cosine(1, 0) = -(direction_cosine(1, 1) = .5 * sin_angle);
direction_cosine(1, 2) = cos_angle;

const auto weight = .5 * edge_length * T;

const mat part_a = shape::stress11(pt_a) * IP(0, 1) * weight;
const mat part_b = shape::stress11(pt_b) * IP(1, 1) * weight;
converter_a = part_a + part_b;
converter_b = .5 * edge_length * (part_a * X1 + part_b * X2);
}

double GCMQ::ResultantConverter::F(const vec& alpha) const { return dot(direction_cosine.row(0), converter_a * alpha); }

double GCMQ::ResultantConverter::V(const vec& alpha) const { return dot(direction_cosine.row(1), converter_a * alpha); }

double GCMQ::ResultantConverter::M(const vec& alpha) const { return dot(direction_cosine.row(0), converter_b * alpha); }

mat GCMQ::form_transformation(const mat& jacobian) {
mat trans_mat(3, 3);

Expand Down Expand Up @@ -142,14 +71,6 @@ int GCMQ::initialize(const shared_ptr<DomainBase>& D) {

access::rw(iso_mapping) = trans(mapping * ele_coor);

const IntegrationPlan edge_plan(1, 2, IntegrationType::GAUSS);
edge.clear();
edge.reserve(4);
edge.emplace_back(ResultantConverter::Edge::A, thickness, ele_coor, edge_plan, iso_mapping);
edge.emplace_back(ResultantConverter::Edge::B, thickness, ele_coor, edge_plan, iso_mapping);
edge.emplace_back(ResultantConverter::Edge::C, thickness, ele_coor, edge_plan, iso_mapping);
edge.emplace_back(ResultantConverter::Edge::D, thickness, ele_coor, edge_plan, iso_mapping);

const IntegrationPlan plan(2, scheme == 'I' ? 2 : 3, scheme == 'I' ? IntegrationType::IRONS : scheme == 'L' ? IntegrationType::LOBATTO : IntegrationType::GAUSS);

const auto diff_coor = form_diff_coor(ele_coor);
Expand Down Expand Up @@ -290,21 +211,15 @@ vector<vec> GCMQ::record(const OutputType P) {
else if(P == OutputType::SP) for(const auto& I : int_pt) data.emplace_back(transform::stress::principal(I.poly_stress * current_alpha));
else if(P == OutputType::SP1) for(const auto& I : int_pt) data.emplace_back(vec{transform::stress::principal(I.poly_stress * current_alpha).at(0)});
else if(P == OutputType::SP2) for(const auto& I : int_pt) data.emplace_back(vec{transform::stress::principal(I.poly_stress * current_alpha).at(1)});
else if(P == OutputType::SINT) data.emplace_back(current_alpha);
else if(P == OutputType::E) for(const auto& I : int_pt) data.emplace_back(I.poly_strain * current_beta);
else if(P == OutputType::E11) for(const auto& I : int_pt) data.emplace_back(I.poly_strain.row(0) * current_beta);
else if(P == OutputType::E22) for(const auto& I : int_pt) data.emplace_back(I.poly_strain.row(1) * current_beta);
else if(P == OutputType::E12) for(const auto& I : int_pt) data.emplace_back(I.poly_strain.row(2) * current_beta);
else if(P == OutputType::EP) for(const auto& I : int_pt) data.emplace_back(transform::strain::principal(I.poly_strain * current_beta));
else if(P == OutputType::EP1) for(const auto& I : int_pt) data.emplace_back(vec{transform::strain::principal(I.poly_strain * current_beta).at(0)});
else if(P == OutputType::EP2) for(const auto& I : int_pt) data.emplace_back(vec{transform::strain::principal(I.poly_strain * current_beta).at(1)});
else if(P == OutputType::EINT) data.emplace_back(current_beta);
else if(P == OutputType::PE) for(const auto& I : int_pt) data.emplace_back(I.poly_strain * current_beta - solve(mat_stiffness, I.poly_stress * current_alpha));
else if(P == OutputType::PE) for(const auto& I : int_pt) { data.emplace_back(I.poly_strain * current_beta - solve(mat_stiffness, I.poly_stress * current_alpha)); }
else if(P == OutputType::PEP) for(const auto& I : int_pt) data.emplace_back(transform::strain::principal(I.poly_strain * current_beta - solve(mat_stiffness, I.poly_stress * current_alpha)));
else if(P == OutputType::RESULTANT) for(const auto& I : edge) data.emplace_back(vec{I.F(current_alpha), I.V(current_alpha), I.M(current_alpha)});
else if(P == OutputType::AXIAL) data.emplace_back(vec{edge[0].F(current_alpha), edge[1].F(current_alpha), edge[2].F(current_alpha), edge[3].F(current_alpha)});
else if(P == OutputType::SHEAR) data.emplace_back(vec{edge[0].V(current_alpha), edge[1].V(current_alpha), edge[2].V(current_alpha), edge[3].V(current_alpha)});
else if(P == OutputType::MOMENT) data.emplace_back(vec{edge[0].M(current_alpha), edge[1].M(current_alpha), edge[2].M(current_alpha), edge[3].M(current_alpha)});
else for(const auto& I : int_pt) append_to(data, I.m_material->record(P));

return data;
Expand Down
18 changes: 0 additions & 18 deletions Element/Membrane/Drilling/GCMQ.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,26 +39,8 @@
class IntegrationPlan;

class GCMQ final : public SGCMQ {
struct ResultantConverter final {
enum class Edge {
A,
B,
C,
D
};

mat converter_a, converter_b;
mat direction_cosine;
ResultantConverter(Edge, double, const mat&, const IntegrationPlan&, const mat&);
[[nodiscard]] double F(const vec&) const;
[[nodiscard]] double V(const vec&) const;
[[nodiscard]] double M(const vec&) const;
};

static constexpr int enhanced_mode = 2;

vector<ResultantConverter> edge;

const mat mat_stiffness, iso_mapping;

mat HT, NT, MT, N, M; // constant matrices
Expand Down
2 changes: 1 addition & 1 deletion Element/Visualisation/vtkBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ class vtkBase {

virtual void Setup();

virtual void GetData(vtkSmartPointer<vtkDoubleArray>&, OutputType = OutputType::U);
virtual void GetData(vtkSmartPointer<vtkDoubleArray>&, OutputType);
virtual mat GetData(OutputType);

virtual void SetDeformation(vtkSmartPointer<vtkPoints>&, double);
Expand Down
Loading

0 comments on commit ca27417

Please sign in to comment.