Skip to content

Commit 150e119

Browse files
elstargazerelstargazer
committed
Merge pull request #15 from elstargazer/fit_ellipse_edited
permit fitting of surfaces not at outer boundary
2 parents 286d4f7 + 2e0cba7 commit 150e119

File tree

3 files changed

+81
-39
lines changed

3 files changed

+81
-39
lines changed

Diff for: ‎config/ConfigurationV2.cfg

+1-1
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
// ===========================================================================
55
// Mesh filename
66
mesh_filename = "meshes/mesh_sph_def_quad_1.inp";
7-
output_folder = "output/2/";
7+
output_folder = "output/temp/";
88
// ___________________________________________________________________________
99
// Body parameters
1010
body_parameters :

Diff for: ‎src/ceres.cc

+27-19
Original file line numberDiff line numberDiff line change
@@ -627,13 +627,10 @@ void StokesProblem<dim>::assemble_system() {
627627
A_Grav_namespace::AnalyticGravity<dim> * aGrav =
628628
new A_Grav_namespace::AnalyticGravity<dim>;
629629
std::vector<double> grav_parameters;
630-
grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep + 0]);
631-
grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep + 0]);
632-
// Note these two core dimensions are fixed to the outer shell and are temporary until a core fitting routine can be re-made
633-
grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep + 0]-system_parameters::depths_rho[0]);
634-
grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep + 0]-system_parameters::depths_rho[0]);
635-
// grav_parameters.push_back(system_parameters::q_axes[1]);
636-
// grav_parameters.push_back(system_parameters::p_axes[1]);
630+
grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep * 2 + 0]);
631+
grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep * 2 + 0]);
632+
grav_parameters.push_back(system_parameters::q_axes[system_parameters::present_timestep * 2 + 1]);
633+
grav_parameters.push_back(system_parameters::p_axes[system_parameters::present_timestep * 2 + 1]);
637634
grav_parameters.push_back(system_parameters::rho[0]);
638635
grav_parameters.push_back(system_parameters::rho[1]);
639636

@@ -1531,12 +1528,17 @@ void StokesProblem<dim>::move_mesh() {
15311528

15321529
// Find ellipsoidal axes for all layers
15331530
std::vector<double> ellipse_axes(0);
1534-
// compute fit to outer boundary
1535-
ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[0]);
1536-
system_parameters::q_axes.push_back(ellipse_axes[0]);
1537-
system_parameters::p_axes.push_back(ellipse_axes[1]);
1538-
std::cout << "a = " << ellipse_axes[0] << " c = " << ellipse_axes[1] << std::endl;
1539-
ellipse_axes.clear();
1531+
// compute fit to boundary 0, 1, 2 ...
1532+
for(unsigned int i = 0; i<system_parameters::sizeof_material_id;i++)
1533+
{
1534+
ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
1535+
system_parameters::q_axes.push_back(ellipse_axes[0]);
1536+
system_parameters::p_axes.push_back(ellipse_axes[1]);
1537+
1538+
std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
1539+
<< " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
1540+
ellipse_axes.clear();
1541+
}
15401542
write_vertices(0);
15411543
}
15421544

@@ -1724,12 +1726,18 @@ void StokesProblem<dim>::setup_initial_mesh() {
17241726

17251727
// Find ellipsoidal axes for all layers
17261728
std::vector<double> ellipse_axes(0);
1727-
// compute fit to outer boundary
1728-
ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[0]);
1729-
system_parameters::q_axes.push_back(ellipse_axes[0]);
1730-
system_parameters::p_axes.push_back(ellipse_axes[1]);
1731-
std::cout << "a = " << ellipse_axes[0] << " c = " << ellipse_axes[1] << std::endl;
1732-
ellipse_axes.clear();
1729+
// compute fit to boundary 0, 1, 2 ...
1730+
std::cout << endl;
1731+
for(unsigned int i = 0; i<system_parameters::sizeof_material_id;i++)
1732+
{
1733+
ellipsoid.compute_fit(ellipse_axes, system_parameters::material_id[i]);
1734+
system_parameters::q_axes.push_back(ellipse_axes[0]);
1735+
system_parameters::p_axes.push_back(ellipse_axes[1]);
1736+
1737+
std::cout << "a_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[0]
1738+
<< " " << " c_"<< system_parameters::material_id[i] <<" = " << ellipse_axes[1] << std::endl;
1739+
ellipse_axes.clear();
1740+
}
17331741
write_vertices(0);
17341742
}
17351743

Diff for: ‎support_code/ellipsoid_fit.h

+53-19
Original file line numberDiff line numberDiff line change
@@ -80,30 +80,64 @@ void ellipsoid_fit<dim>::compute_fit(std::vector<double> &ell, unsigned char bou
8080
std::vector<unsigned int> ind_bnry_col;
8181

8282
// assemble the sensitivity matrix and r.h.s.
83-
for (; cell!=endc; ++cell)
84-
for (unsigned int f=0; f < GeometryInfo<dim>::faces_per_cell; ++f)
85-
{
86-
boundary_ids = cell->face(f)->boundary_indicator();
87-
if(boundary_ids == boundary_that_we_need)
88-
{
89-
for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
90-
if (vertex_touched[cell->face(f)->vertex_index(v)] == false)
83+
double zero_tolerance = 1e-3;
84+
for (; cell != endc; ++cell) {
85+
if (boundary_that_we_need != 0)
86+
cell->set_manifold_id(cell->material_id());
87+
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f) {
88+
if (boundary_that_we_need == 0) //if this is the outer surface, then look for boundary ID 0; otherwise look for material ID change.
9189
{
92-
vertex_touched[cell->face(f)->vertex_index(v)] = true;
93-
for (unsigned int i=0; i<dim; ++i)
94-
{
95-
// sensitivity matrix entry
96-
A(j,i) = pow(cell->face(f)->vertex(v)[i],2);
97-
// r.h.s. entry
98-
b[j] = 1.0;
99-
// if mesh if not full: set the indicator
90+
boundary_ids = cell->face(f)->boundary_indicator();
91+
if (boundary_ids == boundary_that_we_need) {
92+
for (unsigned int v = 0;
93+
v < GeometryInfo<dim>::vertices_per_face; ++v)
94+
if (vertex_touched[cell->face(f)->vertex_index(v)]
95+
== false) {
96+
vertex_touched[cell->face(f)->vertex_index(v)] =
97+
true;
98+
for (unsigned int i = 0; i < dim; ++i) {
99+
// stiffness matrix entry
100+
A(j, i) = pow(cell->face(f)->vertex(v)[i], 2);
101+
// r.h.s. entry
102+
b[j] = 1.0;
103+
// if mesh if not full: set the indicator
104+
}
105+
ind_bnry_row.push_back(j);
106+
j++;
107+
}
108+
}
109+
} else { //find the faces that are at the boundary between materials, get the vertices, and write them into the stiffness matrix
110+
if (cell->neighbor(f) != endc) {
111+
if (cell->material_id() != cell->neighbor(f)->material_id()) //finds face is at internal boundary
112+
{
113+
int high_mat_id = std::max(cell->material_id(),
114+
cell->neighbor(f)->material_id());
115+
if (high_mat_id == boundary_that_we_need) //finds faces at the correct internal boundary
116+
{
117+
for (unsigned int v = 0;
118+
v < GeometryInfo<dim>::vertices_per_face;
119+
++v)
120+
if (vertex_touched[cell->face(f)->vertex_index(
121+
v)] == false) {
122+
vertex_touched[cell->face(f)->vertex_index(
123+
v)] = true;
124+
for (unsigned int i = 0; i < dim; ++i) {
125+
// stiffness matrix entry
126+
A(j, i) = pow(
127+
cell->face(f)->vertex(v)[i], 2);
128+
// r.h.s. entry
129+
b[j] = 1.0;
130+
// if mesh if not full: set the indicator
131+
}
132+
ind_bnry_row.push_back(j);
133+
j++;
134+
}
100135
}
101-
ind_bnry_row.push_back(j);
102-
j++;
103136
}
137+
}
104138
}
105139
}
106-
140+
}
107141
if (ind_bnry_row.size()>0)
108142
{
109143

0 commit comments

Comments
 (0)