classdef Tutorial07_2_PhaseFieldProblemHomogenized < handle

properties (Access = public)

initialGuess

output

end

properties (Access = private)

mesh

boundaryConditions

mat

dissipation

functional

end

methods (Access = public)

function obj = Tutorial07_2_PhaseFieldProblemHomogenized()

obj.init()

obj.defineCase();

obj.createInitialGuess();

obj.createMaterialPhaseField();

obj.createDissipationInterpolation();

obj.createPhaseFieldFunctional()

obj.solvePhaseFieldProblem()

end

end

methods (Access = private)

function init(obj)

close all;

end

function defineCase(obj)

s.mesh.type = '1Elem';

s.bc.type = 'DisplacementTractionY';

s.bc.values = [0:0.001:0.1];

[obj.mesh, obj.boundaryConditions] = BenchmarkManager.create(s);

end

function createInitialGuess(obj)

u = LagrangianFunction.create(obj.mesh,2,'P1');

phi = LagrangianFunction.create(obj.mesh,1,'P1');

%phi = obj.setInitialDamage(phi);

obj.initialGuess.phi = obj.createDamageVariable(phi);

obj.initialGuess.u = u;

end

function phi = setInitialDamage(obj,phi)

isInMiddle = obj.mesh.coord(:,1)>=0.5 & obj.mesh.coord(:,2)==0.5;

fValues = phi.fValues;

fValues(isInMiddle) = 0.01;

phi.setFValues(fValues);

end

function phi = createDamageVariable(obj,phi)

s.type = 'Damage';

s.mesh = phi.mesh;

s.fun = phi;

phi = DesignVariable.create(s);

end

function createPhaseFieldFunctional(obj)

s.energySplit = false;

s.C = obj.mat.C;

s.dC = obj.mat.dC;

s.d2C = obj.mat.d2C;

s.mesh = obj.mesh;

s.dissipation = obj.dissipation;

s.l0 = 0.1;

s.quadOrder = 2;

s.testSpace.u = obj.initialGuess.u;

s.testSpace.phi = obj.initialGuess.phi.fun;

obj.functional = PhaseFieldFunctional(s);

end

function createMaterialPhaseField(obj)

s.fileName = 'CirclePerimeter';

s.mesh = obj.mesh;

s.young = 210;

hm = HomogenizedMaterialsReader(s);

obj.mat.C = @(phi) hm.obtainTensor(phi);

obj.mat.dC = @(phi) hm.obtainTensorDerivative(phi);

obj.mat.d2C = @(phi) hm.obtainTensorSecondDerivative(phi);

end

function createDissipationInterpolation(obj)

s.mesh = obj.mesh;

s.pExp = 2;

Gc = 5e-3;

if s.pExp == 1

cw = 1/2;

elseif s.pExp == 2

cw = 2/3;

end

obj.dissipation.interpolation = PhaseFieldDissipationInterpolator(s);

obj.dissipation.constant = Gc/(4*cw);

end

function solvePhaseFieldProblem(obj)

s.mesh = obj.mesh;

s.initialGuess = obj.initialGuess;

s.boundaryConditions = obj.boundaryConditions;

s.functional = obj.functional;

s.monitoring.set = true;

s.monitoring.type = 'full';

s.monitoring.print = true;

s.tolerance.u = 1e-13;

s.tolerance.phi = 1e-6;

s.tolerance.stag = 1e-6;

s.maxIter.u = 100;

s.maxIter.phi = 300;

s.maxIter.stag = 300;

s.solver.type = 'Gradient';

s.solver.tau = 150;

PFComp = PhaseFieldComputer(s);

obj.output = PFComp.compute();

end

end