From 3db8098ec0695901f2678537aa58fa0bc9d1c1c8 Mon Sep 17 00:00:00 2001 From: Benjamin Spencer Date: Thu, 18 Jun 2015 17:51:05 -0600 Subject: [PATCH] Add PenaltyDirichletBC closes #5268 --- framework/include/bcs/PenaltyDirichletBC.h | 52 ++++++++ framework/src/base/Moose.C | 2 + framework/src/bcs/PenaltyDirichletBC.C | 78 +++++++++++ .../tests/penalty_dirichlet/gold/out.e | Bin 0 -> 44228 bytes .../penalty_dirichlet/penalty_dirichlet.i | 124 ++++++++++++++++++ .../tests/penalty_dirichlet/tests | 7 + 6 files changed, 263 insertions(+) create mode 100644 framework/include/bcs/PenaltyDirichletBC.h create mode 100644 framework/src/bcs/PenaltyDirichletBC.C create mode 100644 modules/solid_mechanics/tests/penalty_dirichlet/gold/out.e create mode 100644 modules/solid_mechanics/tests/penalty_dirichlet/penalty_dirichlet.i create mode 100644 modules/solid_mechanics/tests/penalty_dirichlet/tests diff --git a/framework/include/bcs/PenaltyDirichletBC.h b/framework/include/bcs/PenaltyDirichletBC.h new file mode 100644 index 000000000000..9b154c7bc604 --- /dev/null +++ b/framework/include/bcs/PenaltyDirichletBC.h @@ -0,0 +1,52 @@ +/****************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MOOSE - Multiphysics Object Oriented Simulation Environment */ +/* */ +/* (c) 2010 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/****************************************************************/ + +#ifndef PENALTYDIRICHLETBC_H +#define PENALTYDIRICHLETBC_H + +#include "NodalBC.h" + +class PenaltyDirichletBC; + +template<> +InputParameters validParams(); + +/** + * Boundary condition of a Dirichlet type + * + * Enforce using a penalty method + */ +class PenaltyDirichletBC : public NodalBC +{ +public: + PenaltyDirichletBC(const std::string & name, InputParameters parameters); + +protected: + virtual void computeResidual(NumericVector & residual); + virtual Real computeQpResidual(); + virtual void computeJacobian(); + virtual Real computeQpJacobian(); + virtual void computeOffDiagJacobian(); + + /// The value for this BC + const Real & _value; + + /// The penalty factor + const Real & _penalty; + + // Holds Jacobian entries + DenseMatrix _local_ke; +}; + +#endif /* PENALTYDIRICHLETBC_H */ diff --git a/framework/src/base/Moose.C b/framework/src/base/Moose.C index 8155229e17f4..a2cf50f283be 100644 --- a/framework/src/base/Moose.C +++ b/framework/src/base/Moose.C @@ -58,6 +58,7 @@ // bcs #include "ConvectiveFluxBC.h" #include "DirichletBC.h" +#include "PenaltyDirichletBC.h" #include "PresetBC.h" #include "NeumannBC.h" #include "FunctionDirichletBC.h" @@ -427,6 +428,7 @@ registerObjects(Factory & factory) // bcs registerBoundaryCondition(ConvectiveFluxBC); registerBoundaryCondition(DirichletBC); + registerBoundaryCondition(PenaltyDirichletBC); registerBoundaryCondition(PresetBC); registerBoundaryCondition(NeumannBC); registerBoundaryCondition(FunctionDirichletBC); diff --git a/framework/src/bcs/PenaltyDirichletBC.C b/framework/src/bcs/PenaltyDirichletBC.C new file mode 100644 index 000000000000..988e46c68d01 --- /dev/null +++ b/framework/src/bcs/PenaltyDirichletBC.C @@ -0,0 +1,78 @@ +/****************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MOOSE - Multiphysics Object Oriented Simulation Environment */ +/* */ +/* (c) 2010 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* Under Contract No. DE-AC07-05ID14517 */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/****************************************************************/ + +#include "PenaltyDirichletBC.h" + +template<> +InputParameters validParams() +{ + InputParameters p = validParams(); + p.addRequiredParam("value", "Value of the BC"); + p.addRequiredParam("penalty", "Penalty factor"); + return p; +} + + +PenaltyDirichletBC::PenaltyDirichletBC(const std::string & name, InputParameters parameters) : + NodalBC(name, parameters), + _value(getParam("value")), + _penalty(getParam("penalty")) +{} + +void +PenaltyDirichletBC::computeResidual(NumericVector & residual) +{ + if (_var.isNodalDefined()) + { + dof_id_type & dof_idx = _var.nodalDofIndex(); + _qp = 0; + residual.add(dof_idx, computeQpResidual()); + } +} + +Real +PenaltyDirichletBC::computeQpResidual() +{ + return _penalty * (_u[_qp] - _value); +} + +void +PenaltyDirichletBC::computeJacobian() +{ + if (_var.isNodalDefined()) + { + _qp = 0; + Real jacobian_value = computeQpJacobian(); + dof_id_type dof_index = _var.nodalDofIndex(); + + DenseMatrix & ke = _assembly.jacobianBlock(_var.number(), _var.number()); + _local_ke.resize(ke.m(), ke.n()); + _local_ke.zero(); + + _local_ke(0,0) += computeQpJacobian(); + + ke += _local_ke; + } +} + +Real +PenaltyDirichletBC::computeQpJacobian() +{ + return _penalty; +} + +void +PenaltyDirichletBC::computeOffDiagJacobian() +{ +} diff --git a/modules/solid_mechanics/tests/penalty_dirichlet/gold/out.e b/modules/solid_mechanics/tests/penalty_dirichlet/gold/out.e new file mode 100644 index 0000000000000000000000000000000000000000..4ca9b47de0309b75098d73d4a177a48b69357110 GIT binary patch literal 44228 zcmeHQOOxHkb-uDKiBHR~Bwl3EDJf;UA`WNxuxKW?q!DLGjl_9)IHaX0DFlOyJA)7x zxB_sgnOUT=aQ;EoS*0qKD*pp#k!7l~ae0+}{(!UbCSK$_-DsfkM)!qiGMuiu7wAU! zhtsD|pYA^Wz`c$8UmJ}^^t*)L_vz}pffXl_7aYMi8vO}=7bwj40vGo`#qS%_a1qfJ zzl(`CbCJNE!;?B5Kd@&MpVGe>EN0fko6%!>*G2LN;lzz8MBrSYFxPh}(rfr#5`OE@ ze@an!*j1U%irqv>%fn-D0y+xe-%K&QfMC`#oVQXOYwkvPCx1YBxkBLrC#FE6Lk+u8 z6LzyE>{c3fS-vOWa6Ioy7{Qgp2(BDPaK*4^HnE?Far%gFmI(Yh!Z|O!U>YI|zlHw_ zesAJ;*`9mWGdGI8Fo1XTI)1Nx@GqCfxW3sC{{5*R+6kn=nulJHSYL+G1Ux#UvL&7o zk1pVMaq9Wtzj&ZL<2*#dAK!wI*vW#5O{T-`==H=)d;)_b*?1#dBrB=^yEFU5n~|(Q zmMjS2QF2W2WO|1C%lQ2%ewVl$et8_GVDw(%g>-M=hUkR*1*N~}9_9qWFY`eW$Wa{n z!Cm@e{9Qo&c@)m2To4xVmC}79m35;3>suT6FQ#YkUku{P`z}(%)X}@XCGIc1gX;yv zSKmi*yWS_$lfYbztM8*aE}i;5qWA*htM8+@ecyM0xNkTiB;$m`6HkcW($PH=o=X#4 zMch*!(iZpUh(UZO_LAHXe~IrzyL1%xj__RiXIzPIYJB2bIsVVn`2Ruik#99V<+~yN z<%f87iSn7H5H!7BlV$^9aF>6DzbmBcEK&wDo77}E&-A{u?5=zd_ZOTn2wW#2S)=@s zj8j?2^0>(BdvZEg>&;IeuWxWoA?TNQ?@&4(zxaOA2h!2|=?8a(ATIFuai#;W@<-xN zjeCdUmh@#)+`mw8r5^oe1MUY3Zdcq51y}0gw;FIiQgFNCQdvtr5KpMgQl8>Y@JsS+ z4}QsquDAyjSI9%g-)|FO>Ijo$9PY}e_>*xZAHKsMi9a>&-zafq-oD!q_gIN5>HS_q z+&RLpKC^w)l5j$83<)J1*LAoy)bn4EtQOKy2H?8Vy^A2-6T!-Z@7M&=S`+ta}P`Tqq--Z+J=yP6w}`5ot9*vnD%Dx zw10)PqCTT6ti(QqpLj|#LwS{s@=Rp_ca`Dq3iSu9$aSElWB4h}AK^!GEgi!C>A!^M zDwTtLMpyMmX{l{(dX7g~ImhT;}n$w`*Km!p>)JvwM+iC+>tDp0e%T-A9}3 z?aoAP5o!or8j#%!XEQsPu!k7_vGwR^F>`|?W_VW{DGkrAc?}OtgXeBw`^hN=myzcj z`);x_^czemlj*=H_CyT1*cMc#E<(L?l&gLAr2zcWs zcW?is7hu;5fB(-(4!45&B4PJ2P+wXKzk&>X{?`M8GIr>P&eL3=%Xee8wew{CVLD>( zemg^TAoB-}yvkVIKW3kP_Lto_Sq#ReuI1RyaXlPkW44+_?sLP~Q#T4+-*Q4!X!{80 z$F9nPF-!Zmeb5kxG3r;>W>y4_*$;#Q*^5Hdr`ZoWfEkOs^p!uy*kamADkKE@U?S(}b*=hRw+mc2Lu!824Hu4<`Jbj79!=mlQl*}iv{;u!$a zcm*2K0vp|^j)&`y4|W^!R};)w1CWFsGhIM#7HCG*2eh5Ud*)h)PL0ffO#MPh0Y!zW zg6ac$K`V_*0M^r}4~T{?N@STHl|7dPIWkq2I*Z@ZK>I1q_1bbZ= zSbu9f;h|>-uqy;-ZbV~Y+eiB#Dq&w>#{3Yr@MfhrW;OMzC5W*&^q-X(r<9Gp#wKlp zOwxpnU8@}Qnt;LE&J$q_D*tS0I)2;mY*^^k4yXNRO(qPfrJ+lz4Tnu``L`H#N#ZQ3>_h=6P#u zz0MM4OIu?%o2Y33ey*{J2dXn^0?swI%0X(EJLg9V>h$H(*k@kY~n)hOyuxfW9R&%yG&!Z37;b8gKO-U!w~c&9jCmQrXY1# zadI+HBQ?0)W@^tozkak#a=5*_v%9^uPD`MAIFj9T$^-K#!32{ZMQe-6Xta4|f!CCK zyE~5#)^}*3RS!oLV%affc&gI?^vtLkm06xf;mqQTqiKaQEtKbg`i)c-5JPfKo)Lww zJM9B2hXziV2lJXj8wIQ(#$}*`3Zy{gR#sR0 z&)||A4!SJiCFsG5y$~?TxzPaTND}PUdP2sugUAkI59^of<)jo!tL|y_wQETW7%G}K zS=hdH=uX4P#lCf_d39u7g)J(POI0ru zmi(qqbyKC1`d&%wSmB%|0cx!6sCtifHXm7^JleB1f2Jm_zb^qoMU=LyV2RN+0(-E% z4_N#T*7pzCwY_`RN1yHOZ(0YR?QL=|-C$v88L%=40^r)}8aqCP) zsz>L=M9Bh+VbM@k>WY3T^V5y$phTWgC}kpA{irFPW79q^w54?#M_+xVfK5u)n8GHV z$SGWH+{jwVKC&FS6VFMor6hK2?4v9>{dEqa>(W>O^-CFmBsOBo@z9P-G#czr!4|^C zdP0Kj&E~%6c*!+);k!|VX{7+@Lqt#LNU0L@Qm;M_aSQ<8=2$C*Y0p8clmXbLZ}baA z6iK7L(FY_bI=;-qp)62kW$mEPLluxyTZCK+um~lJ*$hzjl+fk|rQK%wO+Z4!5K9VA zMHux9g~szU^r1a_KDVk>S|1R5cn`yJ4bBPG1$-q(qtt)J!Dh5KIp@5{I2#B)lGQInSd)6Kk2WIInS z?4t=>MUUyLF96PAc%EQ44xO{-MmV-a8La5H%dl{^9dDqk94N~jCR5#Q?*Lr*Ppw%v zS@;WudGsnc6CO%a(8pUr;vPlzk^>*IEd`K^WWhI=gHl&TG@_8?h0eBC7mKCgOZNi-=GKHr%9Fg%Sp*kIitZmR}kfDM-Ru0 z&Z_M1D2-UDD;rTyZI#X1l#ZppI7%TsXSENVL?L#6#4*|xD*Muo5*2&jH!jDmfi@Ei*%W^l#hmhbY9>za$PHNXXt?8G-{Y| zxibMIkd&V4#R%}&#c^)F<(|xGeyqGmUcXIA%7Cim3_@WBD9qbhn8dUyg^@dTaoC$v zETmq*1UUGV`YDdEIEou28#V{+wIW5pjDbWvsjs}A_PKi%C#<(|mcc;dM$Tfk&nk4$LgC>4=D4(ySD}1XPBL3VB>?bd`e~=`|2(Kw)~+ zi#n!Nh|~w<1G8x1B;tzQt`&LUVqprMl1tT|(zH4>0D(v23XE;1bYgAU=~L)nok|pXu-eFMz3 zwR*h@H=`GGqtMlH-RY4c&-$etG!4Bapzhpg(6oKh7_it$PYXA@7|f04@|>Y4H!u-s zt_Tw9YcZdoy7Hj~EOZJN#wayRXk?|>b;H<6>mmY(jc%cwv@Rkr(1)i(Q+PTwg=a%k zcs4ABFTs4uuGktpSf3LHQ#zEFpAn%{xcpXMeL0~YHd!Nqt6x*7_k~X>H1=*WZw z#@$RX^?4`(#j(>R4w@*33=m~A!+I#R_bRG>p`?Jd-B%!E0KT=eb$}Bm*umB|mT+!j zap=xQ<2XFC5!0d;zGMXRuFFw4#|IG;{tZN!Ipw@rRwf1tRse#$W2{}WgCh(~=?guQ zN`0{rfG}%LEkCR-tK5Ig0@oK&f-z3Xrz`|170RYvCr>_rr zmfOA9=FhR7qsY`4WJ%|EPK#x*-p-(UWgoAPyAWQh;@Xvn2rQ2WR>n&TOQTt~4sqNU32$7Luyg;iw64G!rC zBn5@K^C?*=#}z=zJt+rnw&^Yp9)pNT#irhJvqYAk;d3I`wO5p$WlFggP~=)e6e-Jd z8a+1ehU>qpm6*r@$q0l6c@5u&(_QBP3Y)A-^fnXq)oT4R0gN3wPr{pe!XQL+<<1Du zrJc8RmNej+w0bz9LI7G=DOo;T*NDP4`#ER7VN7Aiw2ZMZVQ2)X8qbKSsxzs}_zIbKbe!ZCxgwzMYz3K2i)>TzjI_<*M%{B8 z?70m#v%w0j*`>xvh_j*mxyh8*X%qETF-WTB3={QL0aOR|OAH>r0t@H4K{6E^W19oF zaF&cLq6ceeb9nI$lGwt`Kj~ZC{J2Eb6BbUj+EJ-2QT1M8ld0wsqs&xQGAVJf<^!Lu zW~A_(lGU0_rGzh0D9_Epsc%_Iv#gAAhcQUL$o9k9@kAyz`hDyU`=S00Ye;nEBTXrN zX39Zq#%bJMr=R6iqHBgxpwmPj#_ksj+I5;`Q@t^%dabMHOabW(i<$<&tWcEsb2YQ# z-iX2%-kA2Ho2mx3EWHGKD<07)P}?pof%Tjh@xZFD+fROs<`dMIm6 zBFM9CQYZyLMFfWdR77wXKt%+N0Xq&V=`4pr1frs2E?|K}?r8J%q>!-5Lx3HpcIcW2 zUcF}Y$qW6_TNFKSMuq+wZceJ>{6qa#&-IR~>ug9s$(al*IzFWPaQ(sl`p$#Rrn*rK zCnDsDWLPZxTt1|Ge{26y--83h;uA5V&sKU7(YZ6y*f^`IWqFLQa>#jj{$<{xeD(F^ zb$O5v3JS3~uggU%tckYi_$c;yUG9)@(o<5w=E?zcNK%36tFKqCD^g+c`@a!E?|*Q( z{p-#>Uq#Q5dmZ|%2!mB;=+ Dfq@$A literal 0 HcmV?d00001 diff --git a/modules/solid_mechanics/tests/penalty_dirichlet/penalty_dirichlet.i b/modules/solid_mechanics/tests/penalty_dirichlet/penalty_dirichlet.i new file mode 100644 index 000000000000..36781081591f --- /dev/null +++ b/modules/solid_mechanics/tests/penalty_dirichlet/penalty_dirichlet.i @@ -0,0 +1,124 @@ +[Mesh] + type = GeneratedMesh + dim = 2 + xmin = 0.0 + xmax = 1.0 + ymin = 0.0 + ymax = 1.0 + nx = 1 + ny = 1 + displacements = 'disp_x disp_y' +[] + +[Variables] + [./disp_x] + order = FIRST + family = LAGRANGE + [../] + [./disp_y] + order = FIRST + family = LAGRANGE + [../] +[] + +[AuxVariables] + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[SolidMechanics] + [./solid] + disp_x = disp_x + disp_y = disp_y + [../] +[] + +[AuxKernels] + [./stress_xx] + type = MaterialTensorAux + variable = stress_xx + tensor = stress + index = 0 + [../] + [./stress_yy] + type = MaterialTensorAux + variable = stress_yy + tensor = stress + index = 1 + [../] + [./stress_zz] + type = MaterialTensorAux + variable = stress_zz + tensor = stress + index = 2 + [../] +[] + +[BCs] + [./left_x] + type = PenaltyDirichletBC + variable = disp_x + boundary = left + value = 0.0 + penalty = 1e4 + [../] + [./right_x] + type = PresetBC + variable = disp_x + boundary = right + value = 0.001 + [../] + [./bottom_y] + type = DirichletBC + variable = disp_y + boundary = bottom + value = 0.0 + [../] +[] + +[Materials] + [./solid] + type = Elastic + block = 0 + disp_x = disp_x + disp_y = disp_y + youngs_modulus = 2e4 + poissons_ratio = 0.0 + [../] +[] + +[Executioner] + type = Transient + + #Preconditioned JFNK (default) + solve_type = 'PJFNK' + + petsc_options = '-snes_ksp_ew ' + petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type' + petsc_options_value = '101 hypre boomeramg' + nl_rel_tol = 1e-12 + l_tol = 1e-5 + start_time = 0.0 + dt = 1 + num_steps = 1 +[] + +[Outputs] + file_base = out + output_initial = true + print_linear_residuals = true + print_perf_log = true + [./exodus] + type = Exodus + [../] +[] diff --git a/modules/solid_mechanics/tests/penalty_dirichlet/tests b/modules/solid_mechanics/tests/penalty_dirichlet/tests new file mode 100644 index 000000000000..f215c10afc67 --- /dev/null +++ b/modules/solid_mechanics/tests/penalty_dirichlet/tests @@ -0,0 +1,7 @@ +[Tests] + [./test] + type = 'Exodiff' + input = 'penalty_dirichlet.i' + exodiff = 'out.e' + [../] +[]