Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Prepare elastodynamics #109

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Prepare elastodynamics #109

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 7 additions & 7 deletions apps/contact/src/common.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1903,7 +1903,7 @@ class diffusion_condensed_assembler_nitsche_cells
{
auto fb = make_scalar_monomial_basis(msh, fc, di.face_degree());
auto dirichlet_bf = m_bnd.dirichlet_boundary_func(face_id);
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb, di.face_degree());
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb);
Matrix<T, Dynamic, 1> rhs = make_rhs(msh, fc, fb, dirichlet_bf, di.face_degree());
//ret.block(face_i*fbs, 0, fbs, 1) = mass.llt().solve(rhs);
}
Expand Down Expand Up @@ -2249,7 +2249,7 @@ class diffusion_full_assembler
if (m_bnd.is_dirichlet_face( face_id))
{
auto fb = disk::make_scalar_monomial_basis(msh, fc, di.face_degree());
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb, di.face_degree());
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb);
auto velocity = m_bnd.dirichlet_boundary_func(face_id);
Matrix<T, Dynamic, 1> rhs = make_rhs(msh, fc, fb, velocity, di.face_degree());
svel.block(cbs + i * fbs, 0, fbs, 1) = mass.llt().solve(rhs);
Expand Down Expand Up @@ -2569,7 +2569,7 @@ class diffusion_mix_full_assembler
if (m_bnd.is_dirichlet_face( face_id))
{
auto fb = disk::make_scalar_monomial_basis(msh, fc, di.face_degree());
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb, di.face_degree());
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb);
auto velocity = m_bnd.dirichlet_boundary_func(face_id);
Matrix<T, Dynamic, 1> rhs = make_rhs(msh, fc, fb, velocity, di.face_degree());
svel.block(cbs + i * fbs, 0, fbs, 1) = mass.llt().solve(rhs);
Expand Down Expand Up @@ -2893,7 +2893,7 @@ class contact_full_assembler
if (dirichlet)
{
auto fb = disk::make_scalar_monomial_basis(msh, fc, di.face_degree());
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb, di.face_degree());
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb);
auto velocity = m_bnd.dirichlet_boundary_func(face_id);
Matrix<T, Dynamic, 1> rhs = make_rhs(msh, fc, fb, velocity);//, di.face_degree());
svel.block(cbs + i * fbs, 0, fbs, 1) = mass.llt().solve(rhs);
Expand Down Expand Up @@ -3119,7 +3119,7 @@ class contact_full_assembler_new
auto fb = make_scalar_monomial_basis(msh, face, di.face_degree());
auto dirichlet_fun = m_bnd.dirichlet_boundary_func(face_id);

matrix_type mass = make_mass_matrix(msh, face, fb);// di.face_degree());
matrix_type mass = make_mass_matrix(msh, face, fb);
vector_type rhs = make_rhs(msh, face, fb, dirichlet_fun);// di.face_degree());

sol.block(face_ofs, 0, fbs, 1) = mass.llt().solve(rhs);
Expand Down Expand Up @@ -3448,7 +3448,7 @@ class contact_face_assembler_new
Matrix<T, Dynamic, 1> robin = make_rhs(msh,bfc,fb,bnd.robin_boundary_func(face_id), face_degree);
assert (robin.size() == num_face_dofs);

Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, bfc, fb, face_degree);
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, bfc, fb);

for (size_t i = 0; i < num_face_dofs; i++)
{
Expand Down Expand Up @@ -3518,7 +3518,7 @@ class contact_face_assembler_new
auto fb = make_scalar_monomial_basis(msh, face, di.face_degree());
auto dirichlet_fun = m_bnd.dirichlet_boundary_func(face_id);

Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, face, fb, di.face_degree());
Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, face, fb);
Matrix<T, Dynamic, 1> rhs = make_rhs(msh, face, fb, dirichlet_fun, di.face_degree());

sol.block(face_ofs, 0, fb.size(), 1) = mass.llt().solve(rhs);
Expand Down
4 changes: 2 additions & 2 deletions apps/contact/src/diffusion_nitsche_solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,7 @@ run_hho_diffusion_nitsche_faces(const Mesh& msh,
auto diff = realsol - fullsol;
H1_error += diff.dot(A*diff);

matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);
vector_type u_diff = diff.block(0, 0, cbs, 1);
L2_error += u_diff.dot(mass * u_diff);

Expand Down Expand Up @@ -374,7 +374,7 @@ run_hho_diffusion_nitsche_cells_full(const Mesh& msh,
auto diff = realsol - fullsol;
H1_error += diff.dot(A*diff);

matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);

vector_type u_diff = diff.block(0, 0, cbs, 1);
L2_error += u_diff.dot(mass * u_diff);
Expand Down
16 changes: 8 additions & 8 deletions apps/contact/src/signorini_newton_solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,7 @@ solve_faces(const Mesh& msh, const Function& rhs_fun, const Analytical& sol_fun

H1_increment += du_full.dot(A * du_full);

matrix_type mass = make_mass_matrix(msh, cl, cb);//, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);

vector_type u_diff = du_full.block(0, 0, cbs, 1);
L2_increment += u_diff.dot(mass * u_diff);
Expand Down Expand Up @@ -211,7 +211,7 @@ solve_faces(const Mesh& msh, const Function& rhs_fun, const Analytical& sol_fun
H1_error += diff.dot(Ah*diff);

auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);
vector_type u_diff = diff.block(0, 0, cbs, 1);
L2_error += u_diff.dot(mass * u_diff);

Expand Down Expand Up @@ -419,7 +419,7 @@ solve_faces_hier(const Mesh& msh, const Function& rhs_fun, const Analytical& so

H1_increment += du_full.dot(A * du_full);

matrix_type mass = make_mass_matrix(msh, cl, cb);//, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);

vector_type u_diff = du_full.block(0, 0, cbs, 1);
L2_increment += u_diff.dot(mass * u_diff);
Expand Down Expand Up @@ -473,7 +473,7 @@ solve_faces_hier(const Mesh& msh, const Function& rhs_fun, const Analytical& so
H1_error += diff.dot(Ah*diff);

auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);
vector_type u_diff = diff.block(0, 0, cbs, 1);
L2_error += u_diff.dot(mass * u_diff);

Expand Down Expand Up @@ -643,7 +643,7 @@ solve_cells_full(const Mesh& msh, const Function& rhs_fun, const Analytical& so
diff_sol.block(cell_ofs, 0, num_total_dofs ,1) = du_full;

auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);//, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);
vector_type u_diff = du_full.block(0, 0, cbs, 1);

H1_increment += du_full.dot(Ah * du_full);
Expand Down Expand Up @@ -721,7 +721,7 @@ solve_cells_full(const Mesh& msh, const Function& rhs_fun, const Analytical& so
#endif

auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);

vector_type u_diff = diff.block(0, 0, cbs, 1);
L2_error += u_diff.dot(mass * u_diff);
Expand Down Expand Up @@ -911,7 +911,7 @@ solve_cells_full_hier(const Mesh& msh, const Function& rhs_fun, const Analytica
diff_sol.block(cell_ofs, 0, num_total_dofs ,1) = du_full;

auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);//, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);
vector_type u_diff = du_full.block(0, 0, cbs, 1);

H1_increment += du_full.dot(Ah * du_full);
Expand Down Expand Up @@ -966,7 +966,7 @@ solve_cells_full_hier(const Mesh& msh, const Function& rhs_fun, const Analytica
H1_error += diff.dot(Ah*diff);

auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree());
matrix_type mass = make_mass_matrix(msh, cl, cb);

vector_type u_diff = diff.block(0, 0, cbs, 1);
L2_error += u_diff.dot(mass * u_diff);
Expand Down
2 changes: 1 addition & 1 deletion apps/diffusion/src/diffusion_hho_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ run_hho_diffusion_solver(const Mesh& msh, const size_t degree)
matrix_type ATT = A.block(0,0, cbs, cbs);
vector_type u_cell = fullsol.block(0,0, cbs, 1);
vector_type res = ATT * u_cell - rhs;
matrix_type mm = make_mass_matrix(msh, cl, cb, hdi.cell_degree());
matrix_type mm = make_mass_matrix(msh, cl, cb);
error += res.dot(mm * res);

//auto diff = realsol - fullsol;
Expand Down
8 changes: 4 additions & 4 deletions apps/maxwell/src/maxwell.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -663,7 +663,7 @@ vector_wave_solver(Mesh<T,2,Storage>& msh, size_t order,
auto qps = integrate(msh, cl, 2*chdi.reconstruction_degree());
for (auto& qp : qps)
{
Matrix<T, Dynamic, 2> hphi = rb.eval_curls2(qp.point());
Matrix<T, Dynamic, 2> hphi = rb.eval_curls(qp.point());
Matrix<T, Dynamic, 2> hphi2 = hphi.block(0,0,cb.size(),2);
Matrix<T, Dynamic, 1> ephi = cb.eval_functions(qp.point());
//Matrix<T, Dynamic, 3> rphi = rb.eval_functions(qp.point());
Expand Down Expand Up @@ -1156,7 +1156,7 @@ vector_wave_solver_complex(Mesh<CoordT,3,Storage>& msh, parameter_loader<Mesh<Co
for (size_t i = 0; i < pts.size(); i++)
{
auto phi = cb.eval_functions(pts[i]);
auto cphi = cb.eval_curls2(pts[i]);
auto cphi = cb.eval_curls(pts[i]);
auto ls = phi.transpose()*esolseg;
auto ptid = ptids.at(i);
data_ex[ptid].first += ls(0);
Expand Down Expand Up @@ -1197,7 +1197,7 @@ vector_wave_solver_complex(Mesh<CoordT,3,Storage>& msh, parameter_loader<Mesh<Co
for (auto& qp : qps_f)
{
Matrix<scalar_type, Dynamic, 3> f_phi = fb.eval_functions(qp.point());
Matrix<scalar_type, Dynamic, 3> c_cphi_tmp = cb.eval_curls2(qp.point());
Matrix<scalar_type, Dynamic, 3> c_cphi_tmp = cb.eval_curls(qp.point());
Matrix<scalar_type, Dynamic, 3> c_cphi = disk::vcross(c_cphi_tmp, n);

mass += qp.weight() * f_phi * f_phi.transpose();
Expand All @@ -1217,7 +1217,7 @@ vector_wave_solver_complex(Mesh<CoordT,3,Storage>& msh, parameter_loader<Mesh<Co

Matrix<scalar_type, Dynamic, 3> cphi_tmp = cb.eval_functions(fc_pts[j]);
Matrix<scalar_type, Dynamic, 3> n_x_cphi_x_n = disk::vcross(n, disk::vcross(cphi_tmp, n));
Matrix<scalar_type, Dynamic, 3> ccphi_tmp = cb.eval_curls2(fc_pts[j]);
Matrix<scalar_type, Dynamic, 3> ccphi_tmp = cb.eval_curls(fc_pts[j]);
Matrix<scalar_type, Dynamic, 3> ccphi_x_n = disk::vcross(ccphi_tmp, n);

Matrix<scalar_type, Dynamic, 3> imp = (1./mur)*ccphi_x_n + (jwmu0/Z)*n_x_cphi_x_n;
Expand Down
14 changes: 7 additions & 7 deletions apps/maxwell/src/maxwell_sip_dg.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -222,7 +222,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
for (auto& qp : qps)
{
auto phi = tbasis.eval_functions(qp.point());
auto cphi = tbasis.eval_curls2(qp.point());
auto cphi = tbasis.eval_curls(qp.point());

M += qp.weight() * phi * phi.transpose();
K += qp.weight() * cphi * cphi.transpose();
Expand Down Expand Up @@ -253,7 +253,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
for (auto& fqp : f_qps)
{
auto tphi = tbasis.eval_functions(fqp.point());
auto tcphi = tbasis.eval_curls2(fqp.point());
auto tcphi = tbasis.eval_curls(fqp.point());
auto n_x_tphi = disk::vcross(n, tphi);
auto n_x_tphi_x_n = disk::vcross(n_x_tphi, n);

Expand Down Expand Up @@ -287,7 +287,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
}

auto nphi = nbasis.eval_functions(fqp.point());
auto ncphi = nbasis.eval_curls2(fqp.point());
auto ncphi = nbasis.eval_curls(fqp.point());
auto n_x_nphi = disk::vcross(n, nphi);

Atn += - fqp.weight() * eta_l * n_x_tphi * n_x_nphi.transpose();
Expand Down Expand Up @@ -375,7 +375,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ
{
Matrix<scalar_type, Dynamic, 3> cphi_tmp = cb.eval_functions(fc_pts[j]);
Matrix<scalar_type, Dynamic, 3> n_x_cphi_x_n = disk::vcross(n, disk::vcross(cphi_tmp, n));
Matrix<scalar_type, Dynamic, 3> ccphi_tmp = cb.eval_curls2(fc_pts[j]);
Matrix<scalar_type, Dynamic, 3> ccphi_tmp = cb.eval_curls(fc_pts[j]);
Matrix<scalar_type, Dynamic, 3> ccphi_x_n = disk::vcross(ccphi_tmp, n);

Matrix<scalar_type, Dynamic, 3> imp = (1./mur)*ccphi_x_n - (jwmu0/Z)*n_x_cphi_x_n;
Expand Down Expand Up @@ -513,7 +513,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo
for (auto& qp : qps)
{
auto phi = tbasis.eval_functions(qp.point());
auto cphi = tbasis.eval_curls2(qp.point());
auto cphi = tbasis.eval_curls(qp.point());

M += qp.weight() * phi * phi.transpose();
K += qp.weight() * cphi * cphi.transpose();
Expand All @@ -540,7 +540,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo
for (auto& fqp : f_qps)
{
auto tphi = tbasis.eval_functions(fqp.point());
auto tcphi = tbasis.eval_curls2(fqp.point());
auto tcphi = tbasis.eval_curls(fqp.point());
auto n_x_tphi = disk::vcross(n, tphi);

if (nv.second)
Expand All @@ -558,7 +558,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo
}

auto nphi = nbasis.eval_functions(fqp.point());
auto ncphi = nbasis.eval_curls2(fqp.point());
auto ncphi = nbasis.eval_curls(fqp.point());
auto n_x_nphi = disk::vcross(n, nphi);

Atn += - fqp.weight() * eta_l * n_x_tphi * n_x_nphi.transpose();
Expand Down
26 changes: 16 additions & 10 deletions apps/nonlinear_solid_mechanics/src/contact_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -204,16 +204,19 @@ run_tresca_solver(Mesh<T, 2, Storage>& msh,
return -time * result_type{fx, fy} * exy / (3.0 * (1.0 + material_data.getLambda()));
};

auto solution = [material_data](const disk::point<T, 2>& p) -> result_type
auto solution = [material_data](const disk::point<T, 2>& p, const T time) -> result_type
{
T y = p.y();
T x = p.x();
T exy = std::exp(x * y);
T coeff = 1.0 / (1.0 + material_data.getLambda());

return result_type{x * exy * (1.0 + coeff), y * exy * (-1.0 + coeff)} / 6.0;
return time*result_type{x * exy * (1.0 + coeff), y * exy * (-1.0 + coeff)} / 6.0;
};

auto sol = [solution](const disk::point<T, 2>& p) -> auto
{ return solution(p, 1.0); };

auto s = [rp, material_data](const disk::point<T, 2>& p) -> T
{
T y = p.y();
Expand All @@ -236,7 +239,7 @@ run_tresca_solver(Mesh<T, 2, Storage>& msh,
nl.addBehavior(disk::DeformationMeasure::SMALL_DEF, disk::LawType::ELASTIC);
nl.addMaterialData(material_data);

nl.initial_guess(solution);
nl.initial_guess(sol);

if (nl.verbose())
{
Expand All @@ -254,8 +257,8 @@ run_tresca_solver(Mesh<T, 2, Storage>& msh,
error.h = average_diameter(msh);
error.degree = rp.m_face_degree;
error.nb_dof = nl.numberOfDofs();
error.error_L2 = nl.compute_l2_displacement_error(solution);
error.error_H1 = nl.compute_H1_error(solution);
error.error_L2 = nl.compute_l2_displacement_error(sol);
error.error_H1 = nl.compute_H1_error(sol);

// nl.compute_stress_GP("stress2D_GP_test.msh");
// nl.compute_continuous_displacement("depl2D_cont_test.msh");
Expand Down Expand Up @@ -294,16 +297,19 @@ run_tresca_solver(const Mesh<T, 3, Storage>& msh,
return -time * result_type{fx, 0.0, fz} * exz / (3.0 * (1.0 + material_data.getLambda()));
};

auto solution = [material_data](const disk::point<T, 3>& p) -> result_type
auto solution = [material_data](const disk::point<T, 3>& p, const T& time) -> result_type
{
T z = p.z();
T x = p.x();
T exz = std::exp(x * z);
T coeff = 1.0 / (1.0 + material_data.getLambda());

return result_type{x * exz * (1.0 + coeff), 0.0, z * exz * (-1.0 + coeff)} / 6.0;
return time*result_type{x * exz * (1.0 + coeff), 0.0, z * exz * (-1.0 + coeff)} / 6.0;
};

auto sol = [solution](const disk::point<T, 3>& p) -> auto
{ return solution(p, 1.0); };

auto s = [rp, material_data](const disk::point<T, 3>& p) -> T
{
T x = p.x();
Expand All @@ -322,7 +328,7 @@ run_tresca_solver(const Mesh<T, 3, Storage>& msh,
nl.addBehavior(disk::DeformationMeasure::SMALL_DEF, disk::LawType::ELASTIC);
nl.addMaterialData(material_data);

nl.initial_guess(solution);
nl.initial_guess(sol);

SolverInfo solve_info = nl.compute(load);

Expand All @@ -335,8 +341,8 @@ run_tresca_solver(const Mesh<T, 3, Storage>& msh,
error.h = average_diameter(msh);
error.degree = rp.m_face_degree;
error.nb_dof = nl.numberOfDofs();
error.error_L2 = nl.compute_l2_displacement_error(solution);
error.error_H1 = nl.compute_H1_error(solution);
error.error_L2 = nl.compute_l2_displacement_error(sol);
error.error_H1 = nl.compute_H1_error(sol);

return error;
}
Expand Down
Loading