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.

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

Material model for fracture deformation #1434

Merged
merged 16 commits into from
Sep 28, 2016

Conversation

norihiro-w
Copy link
Collaborator

@norihiro-w norihiro-w commented Sep 23, 2016

This PR adds two models for fracture deformation, linear elastic and Mohr-Coulomb.

Remarks

  • Mohr-Coulomb model is limited to 2D case
  • Fixed-size Eigen matrix/vectors are not considered. let me know if it's really important...

Copy link
Member

@wenqing wenqing left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍

config.checkConfigParameter("type", "LinearElasticIsotropic");
DBUG("Create LinearElasticIsotropic material");

// Youngs modulus
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Young's

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removed

DBUG("Create LinearElasticIsotropic material");

// Youngs modulus
auto& Kn = ProcessLib::findParameter<double>(
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

const?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not necessary. the variable is only passed to MaterialProperties constructor.


DBUG("Use '%s' as Kn parameter.", Kn.name.c_str());

// Poissons ratio
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Poisson's

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

removed

@endJunction
Copy link
Member

@nagelt Could you have a look onto the material models please?

@bilke
Copy link
Member

bilke commented Sep 26, 2016

I guess we are hitting this Qt bug. I also guess the workaround would be to add the proposed #ifndef Q_MOC_RUN as a wrapper around the Boost include in MathTools.h.

@norihiro-w
Copy link
Collaborator Author

@bilke thanks! The workaround solved the error.

@nagelt
Copy link
Member

nagelt commented Sep 27, 2016

Looks good to me. The code can be followed easily.

Fixed size matrices would be sigma/w of (displacement_dim x 1) and the various K/C matrices (displacement_dim x displacement_dim), right? So it's one implementation for MC (2) and two for linear elastic (2,3). @endJunction coming back to Nori's question: How much do we need fixed size when you look at our standard M implementations?

Aside from that I have no further comments on the content.

Copy link
Collaborator

@chleh chleh left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only minor things from my side.

  • Eigen::Ref maybe should be discussed.
  • Regarding fixed-size Eigen matrices we probably should find a common policy at some point. @endJunction only uses fixed-size, if I remember correctly, and here there are only dynamic-size matrices.

@@ -17,7 +17,9 @@
#include <omp.h>
#endif

#ifndef Q_MOC_RUN // to avoid Qt bug, https://bugreports.qt.io/browse/QTBUG-22829
Copy link
Collaborator

@chleh chleh Sep 27, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe add a comment that this is Qt4 specific? ✅
From the link above:

Until you port your code to Qt5 [...]

#ifndef MATERIALLIB_FRACTURE_CREATELINEARELASTICISOTROPIC_H_
#define MATERIALLIB_FRACTURE_CREATELINEARELASTICISOTROPIC_H_

#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
Copy link
Collaborator

@chleh chleh Sep 27, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please move to cpp file. ✅

#ifndef MATERIALLIB_FRACTURE_CREATEMOHRCOULOMB_H_
#define MATERIALLIB_FRACTURE_CREATEMOHRCOULOMB_H_

#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
Copy link
Collaborator

@chleh chleh Sep 27, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also: please move to cpp file. ✅

virtual void computeConstitutiveRelation(
double const t,
ProcessLib::SpatialPosition const& x,
Eigen::Ref<Eigen::VectorXd const> w_prev,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Eigen::Ref wraps arbitrary Eigen expressions. I guess you really intend this. But should we really design our material models s.t. they take Eigen expressions as arguments rather than plain Eigen vectors/matrices?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm using Ref because it can simply take both fixed size and dynamic size Eigen objects. If the performance is important (I doubt because the number of fracture elements are not many), I can change it as Dima did in Solids::MechanicsBase.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds like a cool idea. However, I am not entirely sure if it reduces complexity (fixed & dynamic at the same time) or increases it (now we have to deal with Eigen::Ref, too).
Maybe at some point in the future there should be a discussion about the matrix types used in material models.

const int index_ns = DisplacementDim - 1;
C.setZero();
for (int i=0; i<index_ns; i++)
C(i,i) = _mp.shear_stiffness(t, x)[0];
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Optimization idea: save _mp.shear_stiffness(t, x)[0] to a temporary outside the loop.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i would postpone optimization until we figure out its necessity.

Eigen::Ref<Eigen::VectorXd const> w,
Eigen::Ref<Eigen::VectorXd const> sigma_prev,
Eigen::Ref<Eigen::VectorXd> sigma,
Eigen::Ref<Eigen::MatrixXd> Kep) override;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here the parameter is alled Kep does it have a slightly different meaning now than the C from the base class method?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"e" means elastic, "p" means plastic. Kep means it includes both elastic and plastic parts.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the explanation.

Eigen::Ref<Eigen::VectorXd> sigma,
Eigen::Ref<Eigen::MatrixXd> Kep)
{
if (DisplacementDim == 3)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be more work to static_assert this and to not add a specialization for DisplacementDim == 3?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i am not sure if it works. static_assert means we don't compile fracture stuff at all for 3D, but other fracture model such as elastic one can work with 3D.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see.

Eigen::MatrixXd Ke(2,2);
Ke.setZero();
Ke(0,0) = mat.Ks;
Ke(1,1) = mat.Kn;
Copy link
Collaborator

@chleh chleh Sep 27, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Optimization idea: Maybe there is something possible with Eigen::MatrixXd::Diagonal(...) or similar.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

postponed

}

// check shear yield function (Fs)
double const Fs = std::abs(sigma[0]) + sigma[1] * tan(mat.phi) - mat.c;
Copy link
Collaborator

@chleh chleh Sep 27, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

std::tan, also below. ✅

P const& normal_stiffness, P const& shear_stiffness,
P const& friction_angle, P const& dilatancy_angle,
P const& cohesion)
: normal_stiffness(normal_stiffness), shear_stiffness(shear_stiffness),
Copy link
Member

@endJunction endJunction Sep 28, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Small improvement:
A warning is issued "declaration shadows a field of MaterialLib::Fracture::MohrCoulomb::MaterialProperties" for all constructor's arguments. I usually postfix the arguments with underscore: :white_check_mark:

MaterialProperties(P const& normal_stiffness_) : normal_stiffness(normal_stiffness_) {}


/// converts the given degrees to radians
inline double to_radians(double degrees) {
return degrees*boost::math::constants::pi<double>()/180.;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

FYI: There also is a direct pi/180 constant in boost called radians

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks for the info. but i don't see a reason to change current code.

/// sign function
template <typename T> int sgn(T val) {
return (T(0) < val) - (val < T(0));
}
Copy link
Member

@endJunction endJunction Sep 28, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The sign function is already implemented in boost: #include <boost/math/special_functions/sign.hpp>
I'd prefer an existing function ✅

@norihiro-w norihiro-w merged commit 7373d48 into ufz:master Sep 28, 2016
@norihiro-w norihiro-w deleted the fracture-model branch September 28, 2016 15:09
bilke added a commit to bilke/ogs that referenced this pull request Oct 13, 2016
bilke added a commit to bilke/ogs that referenced this pull request Oct 13, 2016
bilke added a commit to bilke/ogs that referenced this pull request Oct 13, 2016
@ogsbot
Copy link
Member

ogsbot commented Jun 19, 2020

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
7 participants