# Apply integration type constrain condition in FEM

I'm running some simulation use FEM, in my model, I need to apply a constraint condition to the governing equation. My governing equation similar to the diffusion equation as below: $$\frac{\partial c}{\partial t}=\nabla(D\nabla c)\quad\text{on}\quad\Omega$$ and the flux is defined as follow: $$-D\nabla c\cdot n=J$$ And the constraint condition is an integration type condition. For short, I need to constrain the total flux alone the specific surface equals to a constant, which means: $$\int_{\partial\Omega_{1}}(-D\nabla c\cdot n)d\Gamma=\text{Constant}$$ where $$\partial\Omega_{1}$$ is part of my domain's surface, not the total surface!!!

I can write out the weak form for the governing equation(based on N-R method) as follow: $$R_{c}^{I}=\int_{\Omega}\dot{c}N^{I}d\Omega+\int_{\Omega}D\nabla c\nabla N^{I}d\Omega-\int_{\partial\Omega}D\nabla c\cdot n N^{I}d\Gamma$$ and the related stiffness matrix: $$K_{c\dot{c}}^{IJ}=-\frac{\partial R_{c}^{I}}{\partial\dot{c}^{J}}=-\int_{\Omega}N^{J}N^{I}d\Omega$$ $$K_{cc}^{IJ}=-\frac{\partial R_{c}^{I}}{\partial c^{J}}=-\int_{\Omega}D\nabla N^{J}\nabla N^{I}d\Omega+\int_{\partial\Omega_{1}}D\nabla N\cdot n N^{I}d\Omega$$

So how can I apply $$\int_{\partial\Omega_{1}}(-D\nabla c\cdot n)d\Gamma=\text{Constant}$$ to my FEM system?

Someone told me that Lagrange multiplier can solve this problem, but I still don't know how to do it. I'm a newbie of FEM, so some details are greate helpful for me.