dolfin team mailing list archive

["B. Emek Abali" <question151757@xxxxxxxxxxxxxxxxxxxxx>]

New question #151757 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/151757

Dear all, 

I try to use a Newton linearization as a class in python AND implement nonhomogeneous Dirichlet conditions, where am I doing wrong? 

Here is the model problem (cleaned so much that it is not a nonlinear problem anymore but anyhow) the boundary in the first iteration is correct but then it is overwritten with zero (for Newton perturbation this is fine) and I want to put it back in every single iteration. I tried in the def F and def J to do it, but no success! Any better ideas? Or should I copy the whole NonlinearProblem class every time :(

#--------------------------------------------------------------------------
from dolfin import *
parameters["form_compiler"]["cpp_optimize"] = True
class iterate(NonlinearProblem):
    def __init__(self, a, L, bc, exter_B):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
	self.bc = bc
	self.exter = exter_B
    def F(self, b, x):
        assemble(self.L, tensor=b, exterior_facet_domains=self.exter)
	for condition in self.bc: condition.apply(b, x)
    def J(self, A, x):
        assemble(self.a, tensor=A, exterior_facet_domains=self.exter)
	for condition in self.bc: condition.apply(A, x)

mesh = Box(0, 0, 0, 1000.0, 1000.0, 1000.0, 20, 50, 20)
V = VectorFunctionSpace(mesh, 'CG', 1)
D = mesh.topology().dim()
bottom = compile_subdomains('x[1] == 0.0')
top = compile_subdomains('x[1] == 1000.0')
# dummy object into the abstract solver
neumann_domains = MeshFunction("uint", mesh, D-1)
neumann_domains.set_all(0)
shear = Expression('5.0')
bc1 = DirichletBC(V.sub(0), shear, top)
bc2 = DirichletBC(V, (0.0,0.0,0.0), bottom)
bc=[bc1,bc2]
dv = TrialFunction(V)	# increment
w = TestFunction(V)	# test
v = Function(V)		# actual velocity from previous iteration
i, j, k, l = indices(4)
G = 80000.0
lmbd = 120000.0
def vij(v):
	return as_matrix(1.0/2.0*(v[i].dx(j) + v[j].dx(i)), [i,j])

S = as_matrix(lmbd*vij(v)[k,k]*Identity(3)[i,j] + 2.0*G*vij(v)[i,j], [i,j])
form = S[i,j]*w[j].dx(i)*dx 
gain = derivative(form,v,dv)
problem = iterate(gain, form, bc, neumann_domains)
solver = NewtonSolver('gmres') # for nonsymmetric case
solver.parameters['convergence_criterion'] = 'incremental'
solver.parameters['relative_tolerance'] = 1.0e-1
solver.solve(problem, v.vector())
#----------------------------------------

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.